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Abstract 


A  new  "discrete  elements"  (L  )  transport  method  is  derived  and 

N 

compared  to  the  discrete  ordinates  method,  theoretically  and  by 
numerical  experimentation.  The  discrete  elements  method  is  more  accurate 
than  discrete  ordinates  and  strongly  ameliorates  ray  effects  for  the 
practical  problems  studied.  The  discrete  elements  method  is  shown  to  be 
more  cost  effective,  in  terms  of  execution  time  with  comparable  storage 
to  attain  the  same  accuracy,  for  a  one-dimensional  test  case  using 
linear  characteristic  spatial  quadrature.  In  a  two-dimensional  test 
case,  a  vacuum  duct  in  a  shield,  L„  is  more  consistently  convergent 

N 

toward  a  Monte  Carlo  benchmark  solution  than  S  .  using  step 
characteristic  spatial  quadrature.  An  analysis  of  the  interaction  of 
angular  and  spatial  quadrature  in  xy-geometry  indicates  the  desirability 
of  using  linear  characteristic  spatial  quadrature  with  the  L„  method. 

The  discrete  elements  method  is  based  on  discretizing  the  Boltzmann 
equation  over  a  set  of  elements  of  angle.  The  zeroth  and  first  angular 
moments  of  the  directional  flux,  over  each  element,  are  estimated  by 
numerical  quadrature  and  yield  a  flux-weighted  average  streaming 
direction  for  the  element.  (Data  for  this  estimation  are  fluxes  in  fixed 
directions  calculated  as  in  S  .)  The  spatial  quadrature  then  propagates 
the  element  flux  in  this  "steered"  direction.  Since  the  quadrature 
directions  are  not  fixed,  but  are  coupled  to  the  fluxes,  the  method 
strongly  ameliorates  ray  effect.  This  is  verified  using  the  square-in-a- 
square  test  case  originated  by  Lathrop.  A  variety  of  spatial,  angular, 
and  element  quadrature  schemes  are  evaluated  for  both  L„  and  S„.  The 
best  discrete  elements  method  uses  a  hybrid  of  Gauss-Christoffel  polar 
and  composite  3-point  Gauss-Legendre  azimuthal  quadrature. 


I.  Introduction 


The  research  reported  here  has  developed  a  new  numerical  method  for 
solution  of  the  Boltzmann  neutral  particle  transport  equation,  the 
method  of  discrete  elements  (L),  and  tested  its  performance  on  problems 
in  one-  and  two-dimensional  Cartesian  coordinates.  These  problems  have 
been  selected  as  a  proof  of  concept. 

A.  Background 

The  Boltzmann  equation,  in  the  form  of  equation  (1-1),  is  a  balance 
equation  for  the  flux  of  neutral  particles  at  any  point  in  seven¬ 
dimensional  phase  space.  This  flux,  f  ,  is  a  function  of  particle 
position  (x,y,z),  direction  of  particle  motion  ( Q,  4> ) ,  particle  speed  (v) 
or  energy  (E) ,  and  time  (t).  For  the  monoenergetic,  steady  state,  one- 
and  two-dimensional  problems  considered  in  this  research,  the  time  and 
energy  dependences  are  suppressed,  as  are  one  (z)  or  two  (y,z)  of  the 
spatial  dependences. 


+  a  T  ^ 


co  4>+S 

T 


(1-1) 


where 

0  is  the  total  cross-section  for  interaction  (absorption  or  scatter) 
o  is  the  scattering  cross-section 


ft  is  the  direction  of  motion  unit  vector,  ( 6  ,  $  )  in  spherical 
coordinates,  and 

$  is  the  scalar  flux,  and  is  related  to  the  (directionally  dependent) 
flux,  ,  by 

<}>  (x,y,z)  =  /  ^(x,y,z,J!,)  dQ  (1-2) 

4  IT 

Equation  (I— 1)  represents  the  balance  between  loss  rate  (left  side) 
and  gain  rate  (right  side)  which  exists  at  each  point  of  phase  space 
under  steady  state  conditions.  The  first  term  on  the  left  is  the  loss 
rate  of  particles  due  to  divergence,  or  spreading  of  the  flux.  The 
second  term  is  the  loss  rate  due  to  interactions  of  the  particles  with 
the  medium,  either  absorption,  which  destroys  the  particle,  or  scatter, 
which  changes  its  direction  of  motion  (and,  more  generally,  energy, 
although  that  dependence  is  suppressed  here) ,  removing  the  particle  from 
its  original  element  of  phase  space.  The  first  term  on  the  right  is  the 
gain  rate  due  to  particles  which  are  present  at  the  given  space 
location,  but  traveling  in  other  directions,  and  which  scatter,  changing 
to  the  given  direction  of  motion.  The  final  term  is  the  gain  rate  due  to 
creation  of  new  particles  by  any  source  mechanism,  for  example, 
radioactive  decay. 

Actually,  the  (phase  space)  flux,  is  usually  of  interest  only  in 
so  far  as  it  is  needed  to  obtain  its  zeroth  angular  moment,  the  scalar 
flux,  or  its  first  angular  moment,  the  (vector)  current,  J.  The 


scalar  flux  is  of  physical  interest  because  it  determines  rates  of 
reactions,  such  as  fission  or  neutron  activation.  The  current  is  of 
physical  interest  in  determining  the  leakage  rates  of  particles  from  one 
region  to  another,  through  boundary  surfaces. 

Equation  ( I  —  1 )  is  an  integro-dif ferential  equation,  solution  of 
which  is  particularly  difficult.  It  is  solved  (except  for  very  simple 
cases)  by  numerical  means,  such  as  the  method  of  discrete  ordinates, 
rather  than  analytically. 

Since  its  development  in  the  1960's,  the  discrete  ordinates  S., 

N 

method  of  numerical  neutral  particle  transport  has  been  a  mainstay  of 
nuclear  design.  Its  simple,  iterative  computational  structure  supports 
computer  solution  of  problems  with  large  spatial  grids  and  multiple 
energy  groups  without  prohibitive  storage  requirements.  This  is 
accomplished  by  treating  the  spatial  and  angular  integrations,  or 
quadratures,  independently.  Finite  element  methods  which  couple  the 
spatial  and  angular  representations  of  the  flux  have  the  potential  for 
accurate  solutions  with  coarse  meshes,  but  are  complicated  and  very 
costly  in  storage. 

Discrete  ordinates  methods  have  serious  drawbacks,  however.  The  use 
of  discretized  angular  and  spatial  representations  inevitably  entails 
truncation  errors.  These  often  take  the  form  of  random  errors  which 
limit  the  accuracy  of  the  method,  but  when  the  truncation  errors  appear 
systematically,  they  can  cause  results  which  are  qualitatively 
incorrect.  Such  systematic  errors  have  been  known  as  ray  effects,  since 
the  discrete  ordinates  method  represents  particle  motion  only  in  fixed 
directions,  or  "rays".  Another  limitation  of  discrete  ordinates  is  that 
it  provides  poor  accuracy  for  some  types  of  problems  which  are  important 


in  nuclear  design.  Streaming  ducts  in  shields  cause  such  difficulties 
because  the  discrete  ordinates  angular  representation  is  not  well 
adapted  to  a  distribution  of  flux  which  is  strongly  peaked  along  the 
axis  of  the  duct.  If  the  duct  is  not  aligned  along  one  of  the  quadrature 
set  directions,  the  discrete  ordinates  calculation  may  be  very 
insensitive  to  the  presence  of  the  duct.  Monte  Carlo  methods  can  be  used 
in  such  cases,  but  are  computationally  costly.  The  performance  of  the  Sv, 
method  is  sensitive  to  the  choice  of  the  angular  quadrature  direction 
set,  in  a  problem-dependent  way.  One  of  the  limitations  of  the  method  is 
that  when  different  quadrature  sets  give  widely  varying  results,  it  is 
difficult  to  have  much  confidence  in  any  of  them. 

The  discretized  spatial  representation  can  also  lead  to  systematic 
errors  which  result  in  qualitatively  incorrect  solutions,  here  called 
"quasi-ray  effect".  Recent  progress  has  been  made  in  developing  higher 
order  spatial  quadrature  schemes  which  reduce  these  truncation  errors. 

Various  schemes  have  been  proposed  for  improving  the  performance  of 
discrete  ordinates  calculations.  Some  of  these  are  reviewed  in  the  body 
of  this  report. 

B.  Statement  of  the  Problem 

The  objective  of  this  research  is  to  derive  the  discrete  elements 
equations,  and  to  develop,  implement,  and  evaluate  the  performance  of 
the  discrete  elements  method.  The  method  should  use  an  improved  angular 
representation  which  couples  angular  and  spatial  quadrature  in  order  to 
reduce  truncation  errors  and  ray  effects,  while  retaining  the 

computational  simplicity  of  the  S„  method. 
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C.  Scope 

The  scope  of  the  research  includes  development  and  demonstration  of 
the  discrete  elements  method  and  analysis  of  its  computational  costs  in 
terms  of  computer  storage  and  execution  time.  Its  performance  is 
compared  to  discrete  ordinates  calculations  with  various  quadrature 
sets.  The  method  is  tested  with  various  angular  and  spatial  quadratures 
and  the  optimum  scheme  identified.  Two  1-dimensional  and  two 
2-dimensional  test  problems  are  used.  These  test  problems  incorporate 
the  following  assumptions: 

-  Geometries:  1-D  and  2-D  Cartesian 

-  Sources: 

Uniformly  distributed  (by  region) 

Isotropic 

-  Media: 

Uniform  by  region 

Isotropic 

Non-Multiplying 

-  Scatter:  Isotropic  (in  Lab  Frame) 

-  Energy  Dependence:  Monoenergetic  (i.e.  One  Group) 

-  Time  Dependence:  Steady  State 

Since  the  discrete  element  method  retains  much  of  the  computational 
structure  of  discrete  ordinates  methods,  the  extension  to  multiple 
energy  groups,  anisotropic  scatter,  and  time  dependent  problems  is 
immediate  and  so  need  not  be  explicitly  demonstrated  here. 


D.  General  Approach  and  Sequence  of  Presentation 

The  discrete  elements  equations  are  derived,  in  chapter  II,  by 
discretizing  the  integro-dif ferential  form  of  the  Boltzmann  transport 
equation  over  elements  of  solid  angle.  The  computer  algorithm  for  the 
method  is  then  developed  as  an  extension  of  the  discrete  ordinates 
algorithm.  Angular  quadrature  sets  for  S„  methods  are  reviewed  in 
chapter  III,  and  corresponding  quadrature  grids  for  L„  are  proposed. 
Three  quadrature  rules  are  considered  in  chapter  IV.  These  are  used 
within  each  element  of  angle  for  coupling  the  angular  and  spatial 
quadratures.  Potential  advantages  and  disadvantages  of  the  three  methods 
are  evaluated.  Chapter  V  reviews  the  spatial  quadrature  schemes  that  are 
to  be  used  in  the  testing  in  later  chapters.  The  review  is  conducted  in 
detail  in  order  to  support  the  later  analysis  of  the  performance  of  the 
discrete  elements  method  using  these  schemes.  Smoothness,  pointwise 
accuracy,  global  accuracy,  and  tendencies  toward  systematic  error  are 
considered.  Having  established  the  angular  and  spatial  quadratures 
available,  the  interaction  of  the  two  forms  of  quadrature  is 
investigated  in  chapter  VI.  Formal  definitions  are  given  for  ray  effect 
and  quasi-ray  effect,  and  examples  of  each  are  given.  In  evaluating  the 
performance  of  the  LN  method,  it  is  important  to  be  able  to  distinguish 
between  faults  of  the  angular  representation  and  quadrature  and  faults 
of  the  spatial  discretization  and  quadrature. 

Since  there  are  many  possible  combinations  of  quadrature  schemes  to 
be  employed  as  a  discrete  elements  method,  chapter  VII  defines  a 
notational  system  for  identifying  the  different  schemes,  and  discusses 
the  possibilities  for  optimization  of  the  method.  Chapter  VIII  then 
analyses  the  computational  costs  of  implementing  these  various  schemes. 
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The  discrete  elements  method  was  programmed  and  tested  in  both  one- 
and  two-dimensional  Cartesian  geometry.  Chapter  IX  presents  an  analysis 
of  the  results  for  the  one-dimensional  problems:  a  pure  absorber  (no 
scatter,  no  source)  with  an  isotropic  flux  incident  on  one  side;  and  a 
highly  scattering  source  region  with  a  low  scatter  shield  and  vacuum 
outer  boundary.  Chapter  X  treats  two  2-dimensional  problems:  a  square 
source  in  a  larger  square  absorber;  and  a  shielded  source  with  a  vacuum 
duct  in  the  shield.  Both  chapters  present  the  results  of  calculations 
with  various  L  schemes,  as  well  as  with  competing  S„  schemes,  and, 
in  most  cases,  analytic  or  Monte  Carlo  benchmark  solutions.  The  methods 
are  analyzed  for  accuracy,  convergence,  and  cost-effectiveness.  Causes 
of  the  observed  errors  are  also  considered.  Chapter  XI  summarizes  the 
conclusions  drawn  in  the  previous  chapters  and  presents  recommendations 
for  use  of  the  discrete  element  method  and  for  future  research. 
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II.  Derivation  of  the  Discrete  Elements  Transport  Method 


This  chapter  presents  the  derivation  of  a  new  method  of 
discretizing  and  solving  the  integro-dif ferential  form  of  the  Boltzmann 
neutral  particle  transport  equation.  The  method  has  some  features  in 
common  with  the  discrete  ordinates  S„  method,  but  differs  from  it  in  two 
major  ways.  First,  the  Boltzmann  equation  is  discretized  over  a 
collection  of  elements  of  solid  angle,  rather  than  into  a  point-set  of 
fixed  directions  as  in  SN*  Second,  the  method  of  approximating  and 
solving  the  discrete  element  equations  explicitly  couples  the  angular 
dependence  of  the  directional  flux  into  the  spatial  quadrature  scheme, 

4 

so  that  the  directions  of  streaming  are  interpolated  to  better  represent 
the  motion  of  the  particles  in  each  angular  element.  A  motivation  for 
this  scheme  is  the  anticipation  that  the  steered  streaming  of  the 
discrete  elements  method  should  ameliorate  the  ray  effects  seen  in  the 
discrete  ordinates  method  as  a  result  of  its  fixed  streaming  directions. 

A.  Angular  Coordinates  and  Elements  of  Angle 

Three  related  angular  coordinate  systems  are  used  here.  These  are: 

1  -  polar  coordinates  (0,$)  where  0  is  the  polar  angle, 
measured  from  the  z-axis,  and  $  is  the  azimuthal  angle,  measured  from 
the  x-axis  toward  the  y-axis. 

2  -  polar  cosine  /  azimuthal  angle  coordinates  (t,$)  where 
x  ■  cos  (0)  and  $  is  the  azimuthal  angle,  as  above. 


3  -  direction  cosine  coordinates  (y,n)  where  U  is  the 


x-direction  cosine  and  n  is  the  y-direction  cosine. 
These  coordinates  are  related  by 


T 

=  cos 

(0) 

(II-l) 

V 

=  sin 

(0) 

cos  {$) 

-Vl-X2 

cos  (<j>) 

(II-2) 

n 

■  sin 

(0) 

sin  ($) 

■  V1  ■  t  2 

sin  ($) 

(II-3) 

The  transport  equation  will  be  discretized  into  a  system  of 
equations,  each  describing  the  flux  over  an  element  of  angle.  These 
elements  of  solid  angle  are  like  wedges  or  cones  which  together  form  the 
unit  sphere  of  solid  angle.  The  angular  mesh  used  here  is  similar  to  a 
product  quadrature  set,  and  is  formed  by  dividing  the  polar  and 
azimuthal  coordinates,  independently,  into  segments.  This  is  perhaps 
most  easily  visualized  as  a  partitioning  of  the  surface  of  a  globe  along 
lines  of  latitude  and  of  longitude.  This  partition  maps  as  a  rectangular 
mesh  on  the  (x,$)  plane,  where  the  weight  of  an  element  is  proportional 
to  its  area.  The  element  weight  can  also  be  expressed  as  the  product  of 
the  x-weight  and  the  $ -weight; 

W  -  /  sin  (9)  d9  d$  /  2ir  =  At  A*  /  2  x  (II-4) 

K  *  1  ^  K  1 

m 

where  D  is  the  domain  of  angle  of  the  m'th  discrete  element  and  where 
m 

the  normalization  factor  is  2  x  rather  than  4x,  since  only  the  upper 

hemisphere  of  directions  need  be  considered  (due  to  the  z-symmetry  of 

xy-geometry) .  Except  where  explicitly  needed,  the  indices,  k  and  1,  will 

» 

be  expressed  as  a  single  index,  m,  ranging  from  1  to  M. 
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B.  Discretization  of  the  (Angular)  Flux 


In  order  to  discretize  the  the  transport  equation,  it  is  useful  to 

first  treat  the  flux,  (x,y,x,$).  The  flux  may  be  considered  as  a  union 

of  M  functions,  ij>  (x,$),  each  over  a  domain  restricted  to  a  single 

m 

angular  mesh  cell,  D  ,  of  area  W  .  One  reason  for  doing  this  is  that 

m  m 

these  functions  could  be  approximated  by  a  truncated  power  series 

expansion.  This  representation  of  the  flux  as  a  piecewise  polynomial  is 

similar  to  a  finite  elements  method  approach.  Differences  are  (1)  finite 

elements  typically  uses  basis  elements  which  overlap  in  support  and  (2) 

finite  elements  typically  uses  differentiability  constraints,  resulting 

in  coupling  among  the  discretized  equations,  and  thus  requiring  the 

assembly  of  global  matrices  and  solution  of  linear  algebra  problems.  The 

discrete  elements  method,  however,  uses  basis  elements  that  are  discrete 

in  domain  (non-overlapping)  and  does  not  require  differentiability,  or 

even  continuity,  at  the  boundaries  between  elements,  thus  avoiding  the 

complications  and  computational  difficulties  of  finite  elements. 

Each  flux,  ,  can  be  characterized  by  its  angular  moments.  The 
m 

zeroth  and  first  moments  are  central  to  the  discrete  element  method. 
Notation  for  these  flux  moments  is  defined  as  follows: 


F  =  (1/w  )  f  dQ 

m  m  D 

m 


J  =  (1/WJ  /  n  1|>  (x,<|>)  dfi 

— m  m  _ 

D 

m 


where  the  directional  unit  vector,  Q,  is  defined  as 


yi+nj+xk 
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F  may  be  physically  interpreted  as  the  average  flux  of  particles 
m 

traveling  in  directions  which  fall  within  the  m'th  element  of  solid 
angle.  Similarly,  the  vector  quantity,  n$(ft),_is  a  (directionally 
dependent)  current.  It  bears  the  same  relationship  to  the  (scalar) 
current  vector,  J,  of  diffusion  theory  as  does  the  (directional)  flux  of 
transport  theory  to  the  (scalar)  flux  of  diffusion  theory.  Thus,  is 
the  average  current  vector  for  the  m'th  element  of  solid  angle. 


C.  Discretization  of  the  Transport  Equation 


The  transport  equation  is  usually  written  in  the  form 


■  co*+S 
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where  c  *  o  /a  and  ♦  is  the  scalar  flux.  Since  the  current  is 
s  t 

physically  meaningful,  as  discussed  in  the  previous  section,  the  first 
step  in  the  discretization  is  to  bring  the  unit  direction  vector  back 
inside  the  divergence  operator: 


7»  (fl'P)  +  ■  ccf+S 


( 1 1—9) 


In  order  to  discretize  this  transport  equation,  it  is  averaged  over 

each  element  of  solid  angle.  More  precisely,  the  zeroth  angular  moment 

of  the  equation  is  taken,  over  the  domain  D  ,  and  normalized  by  the 

m 

weight,  H  : 

tn 


w  fL  7*  cfi'PJdn  +  77-  /n  o'Mn  -  cx*+s 


(11-10) 


Since  the  spatial  divergence  and  angular  moment  operators  apply  to 
independent  variables,  the  divergence  operator  may  be  taken  outside  the 


integral.  Recalling  the  definitions  of  F  and  J  ,  the  result  is 

m  — m 


V  •  J  +  o  F  =  c  o  $  +  S 
— m  m 


(II-ll) 


Two  simplified  forms  of  this  equation  are  convenient.  In  a  vacuum,  it 


reduces  to 


7*  j  =  s 

-m 


(11-12) 


and  otherwise,  since  the  cross-section  is  non-zero,  it  may  be  written 


-  V*  J  +  F  =  Q 
0  —to  m  x 


(II-13) 


where  Q  is  the  effective  source  term,  with  dimensions  of  flux,  defined 


Q  -  c*  +  - 
o 


(11-14) 


The  scalar  flux  is  eliminated  from  this  system  of  equations  by 


expressing  it  in  terms  of  the  F 


♦  -  /  *  d  n 

4ir 


(11-15) 


Q-T+c^WF 
*  d  mm 


(11-16) 


Equations  (11-16)  and  (11-13),  or  (11-12)  in  vacuum,  are  the 
discrete  elements  equations.  These  equations  are  exact,  which  is  to  say 
that  no  approximations  have  been  made  yet.  If  they  were  solved 
analytically,  the  exact  scalar  flux  and  current  would  be  obtained, 
although  the  directional  flux,  ’J'  ,  would  not  be  fully  specified. 
Approximations  to  the  second  or  higher  angular  flux  moments  could  be 
devised,  but  the  first  moment  is  determined  exactly  by 


i  -  L  \ 


(11-17) 


Although  equation  (11-16)  assumes  isotropic  scatter,  equation  (11-17) 
could  be  used  to  include  linearly  anisotropic  scatter  in  a 
straightforward  manner.  For  simplicity  in  implementing  and  evaluating 
the  discrete  elements  method,  isotropic  scatter  is  assumed. 


D.  The  Flux-Weighted  Mean  Angle 


The  discrete  elements  equations  appear  very  much  like  the  discrete 

ordinates  equations.  In  fact,  if  an  assumption  were  made  that  the 

discrete  fluxes  and  currents  were  related  by  J  *  ft  F  where  the 

— m  m  m 

direction  (or  "ordinate"),  n  ,  is  a  fixed  direction  established  by  the 

m 

choice  of  angular  quadrature  set,  then  the  discrete  elements  equations 
would  reduce  to  the  discrete  ordinates  equations.  Without  this 


assumption,  the  relationship  of  and  F^  defines  the  element  mean 
streaming  direction  to  be 


J_  /  F_ 


(11-18) 


Substituting  equations  (I I— 5)  and  (II-6),  and  canceling  the  weights,  W  , 

m 


n(ft)  dfi 


/  <Mft)  dft 


(11-19) 


The  components  of  ft  are  designated  by  0  ,  <j>  ,  y  ,  n  ,  and  t  . 

m  m  m  m  m  m 

A 

Equation  (11-19)  shows  ft  to  be  the  element  flux-weighted  mean  angle.  If 

m 

these  angles  were  somehow  known,  then  the  well-developed  techniques  for 
solving  the  discrete  ordinates  equations  could  be  used  to  solve  the 

X 

discrete  elements  equations.  Numerical  methods  of  approximating  the  ft 

m 

are  considered  next. 


atinq  the  Flux-Weighted  Streaming  Directions 


The  flux-weighted  mean  angles,  ft  ,  are  approximated  by  using 

m 

numerical  quadrature  rules  to  approximate  the  integrals  in 

equation  (11-19).  Separation  of  variables  in  ( t  ,  «J>  )  coordinates  is 

assumed  within  each  discrete  element,  and  the  integrals  in  the  two 

coordinates  are  done  independently.  Thus,  for  each  element,  there  is  a 

flux-weighted  mean  t  ,  and  a  flux-weighted  mean  <f>  : 

m  m 


/  xty( t )  dx 
D 

m _ 

/  <Mt)  dt 


(11-20) 


f  <fiip(<P)  d4 1 


f  (<f>)  d$ 


(11-21) 


These  one-dimensional  integrals  are  then  approximated  by  a 


quadrature  rule,  e.g.  Simpson's  rule.  The  application  of  Simpson's  rule 
requires  knowledge  of  the  flux  in  three  fixed  directions  distributed  in 
azimuth  and  three  fixed  directions  distributed  polarly,  within  each 
element.  Figure  II-l  shows  these  fixed  "auxiliary"  directions,  labelled 
as  north,  south,  east,  west,  and  center.  This  nomenclature  is  chosen  to 
be  analogous  to  the  sides  of  a  rectangle  of  latitude  and  longitude  on 
the  surface  of  a  globe,  eastward  being  increasing  <J>,  and  northward  being 
increasing  t. 


t  T 


Fig.  II-l:  Discrete  Element  Auxiliary  Directions  for  Simpson's  Rule 

These  directions  are  fixed  by  the  choice  of  the  discrete  element 
angular  mesh  and  the  choice  of  Simpson's  rule.  The  flux  in  these  fixed 
directions  therefore  obeys  the  discrete  ordinates  transport  equations 
and  may  be  calculated  using  conventional  SN  spatial  quadrature.  The 
source  term  for  these  auxiliary  calculations  is  the  same  Q  as  is  used  in 
the  discrete  elements  equations.  The  algorithm  for  this  method  is 
considered  in  the  next  section. 

First,  an  alternative  view  of  this  approximation  scheme  is 
considered.  If  the  flux  is  expanded  in  a  Taylor's  series  (in  angle) 
about  the  element  center,  it  is  expressed  as  a  series  in  powers 
of  (t-t  )  and  Retaintng  terms  through  quadratic,  including  the 


North 

- e - 

( i  o  Center 

- o 

South 

— > 


bilinear  term,  and  performing  the  integrations  in  equation  (11-19) 
analytically,  it  may  be  shown  that  the  bilinear  term  vanishes  from  both 
integrals.  (Thus,  the  assumption  of  separability  is  actually  required 
only  if  terms  higher  than  quadratic  are  retained.)  The  use  of  Simpson's 
rule  is  equivalent  to  finding  a  quadratic  flux  which  equals  the 
calculated  fluxes  at  the  five  auxiliary  directions  and  exactly 
evaluating  the  flux-weighted  mean  streaming  direction  for  that  flux. 

F.  The  Discrete  Elements  Method  Algorithm 

The  discrete  elements  algorithm  is  based  as  closely  as  possible 

upon  that  used  in  discrete  ordinates.  The  essential  feature  of  S„  is  the 

separation  of  the  angular  and  spatial  quadratures  in  an  iterative 

scheme.  The  source  term,  Q,  is  assumed  known  (from  the  previous 

iteration,  or  by  a  guess),  and  the  spatial  quadrature  is  performed  as  a 

walk  through  a  mesh  of  space  cells,  using  the  fixed  directions  of  the 

angular  quadrature  set.  As  the  directional  fluxes  are  found,  they  are 

folded  into  a  new  scalar  flux  at  each  point,  thus  performing  the  angular 

quadrature  to  get  an  improved  scalar  flux.  After  the  spatial  walk  is 

completed,  the  source  term,  Q,  is  updated  in  preparation  for  the  next 

iteration.  This  scheme  minimizes  storage  requirements  so  that  real 

engineering  problems  can  be  fit  into  the  computer.  The  discrete  elements 

algorithm  retains  this  structure,  with  its  computational  advantages,  but 

adds  the  steps  necessary  to  approximate  the  (no  longer  fixed)  streaming 

directions  and  use  them  in  the  spatial  quadrature  of  the  element  fluxes, 

F  .  Once  the  streaming  direction  is  found,  the  spatial  quadrature  of  F 
m  m 

uses  the  same  methods  as  in  S^. 
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The  discrete  elements  iterative  scheme  is: 

1.  Perform  space  quadrature  in  fixed  directions  to  get 
auxiliary  directional  fluxes. 

2.  Perform  angular  quadrature  within  each  element  of  angle  to 

A 

obtain  the  steered  ft 

m 

3.  Perform  space  quadrature  in  the  (adjusted)  directions  to 

get  discrete  element  ("main")  directional  fluxes. 

4.  Accumulate  the  (weighted)  main  fluxes  to  get  the  scalar 
flux  at  each  space  cell. 

; 

5.  Update  the  source  term,  Q,  and  return  for  next  iteration, 

if  the  results  are  not  yet  converged. 


The  actual  algorithm  combines  steps  1  through  4  in  order  to 
minimize  storage  requirements,  as  follows: 


Repeat 

For  each  angular  element 
For  each  space  cell 

1.  Step  the  associated  auxiliary  fluxes 

across  the  cell 

2.  Find  the  mean  direction  for  the  main  flux 

3.  Step  the  main  flux  across  the  cell 

4.  Accumulate  the  contribution  to  the  scalar 

flux  in  the  "new-flux"  array 
Next  space  cell 
Next  angular  element 

For  each  space  cell 

1.  Test  for  convergence 

2.  Update  the  "old"  scalar  fluxes  to  the  "new" 

3.  Re-evaluate  the  source,  Q,  for  the  cell 
Next  space  cell 


Until  convergence  criteria  are  met. 


G.  Summary 


The  discrete  elements  equations  and  algorithm  have  been  developed 

in  this  chapter.  The  method  is  similar  to  a  finite  elements  method  in 

that  the  flux  is  represented  by  piecewise  polynomial  functions.  The 

advantages  of  finite  elements  may  be  expected,  specifically,  better 

performance  for  problems  with  strong  local  sources  and  absorbers,  since 

the  resulting  ill-behaved  flux  functions  are  better  approximated  by 

the  piecewise  basis.  Unlike  the  finite  elements  method,  however,  the 

elements  are  discrete  and  are  solved  by  an  iterative  spatial  and  angular 

mesh  walk  rather  than  simultaneously  with  matrix  algebra.  In  this 

respect,  the  method  retains  the  practicality  of  the  discrete  ordinates 

method.  Unlike  S„,  however,  the  spatial  and  angular  quadratures  are 
N 

explicitly  coupled  through  the  flux-weighted  mean  streaming  angles.  In 

effect,  this  coupling  steers  the  streaming  in  the  average  directions  of 

particle  flow  and  should  strongly  ameliorate  the  ray  effects  caused  by 

using  fixed  streaming  directions  in  S  . 

N 

Before  the  L *  method  can  be  implemented  and  evaluated  by  comparison 
N 

with  S„,  the  details  of  the  angular  and  spatial  quadrature  methods  must 

be  filled  in.  Chapter  III  covers  angular  quadrature:  choice  of  streaming 

directions  and  weights  in  S  .  and  choice  of  angular  mesh  in  L„.  Chapter 

N  N 

IV  covers  the  angular  quadrature  used  within  each  element  to  evaluate 
the  mean  streaming  direction.  Chapter  V  presents  the  spatial  quadrature 
methods  used  to  step  the  fluxes  across  each  spatial  grid  cell.  Chapter 
VI  considers  the  interaction  of  the  spatial  and  angular  discretized 
quadratures  and  the  inherent  deficiencies  of  discretization,  such  as  ray 
effect. 
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III.  Angular  Quadrature  in  the  and  Methods 


In  the  discrete  ordinates  method,  a  quadrature  set  is  a  set  of 

A 

directions,  and  weights,  W^,  used  to  approximate  angular  integration 

by  the  following  quadrature  rule: 


f  f(fi)  dft 

D 

m 


)  w  f  (n  ) 

m  m 
m 


(III-l) 


where  f(ft)  is  any  function,  for  example,  .  If  equation  (III-l)  is  exact 

for  f(&)  =  pn  ,  then  the  quadrature  set  is  said  to  match  the  n'th 

angular  moment  (with  respect  to  the  x-axis,  in  this  case).  This  chapter 

reviews  some  of  the  quadrature  sets  proposed  or  in  use  for  S  ,  and 

N 

defines  quadrature  sets  for  use  in  one-  and  two-dimensional  (slab  and 

xy)  geometry  with  the  L„  method. 

N 


A.  Global  Basis  Functions  and  Gaussian  Quadrature 

SN  conventionally  uses  quadrature  sets  which  match  as  many  angular 
moments  of  the  integrand  as  possible.  (Various  symmetry  constraints  are 
usually  imposed,  as  well.)  If  the  integrand  is  well-behaved,  in  that  it 
has  high  orders  of  differentiability,  then  it  can  be  expanded  in  a  power 
series  in  the. direction  cosines.  The  quadrature  set  then  integrates, 
exactly  (or  nearly  exactly),  as  many  of  the  terms  as  possible.  Gaussian 
quadrature  is  an  example  in  one  dimension.  This  works  well  if  the  flux 
is  well  represented  by  a  truncated  power  series  in  the  direction 


cosines.  An  alternative  approach  is  discussed  next. 


B.  Local  Basis  Functions  and  Composite  Quadrature 


models  the  flux  with  piecewise  low  order  polynomials  in  one  or 
two  angular  coordinates,  x  and/or  ^ .  These  are  integrated  by  low-order 
methods  such  as  Simpson's  rule.  The  summation  of  the  integrals  over  the 
pieces  constitutes  a  composite  quadrature.  If  the  integrand  of  a 
composite  integration  is  actually  well-behaved,  the  composite 
integration  works  well,  although  a  Gaussian  would  be  more  cost 
effective.  But,  if  the  integrand  is  ill-behaved  (in  that  its  low-order 
derivatives  are  small  but  its  high-order  derivatives  are  large  or 
unbounded),  the  composite  rule  improves  in  accuracy  as  more  subintervals 
are  used,  while  Gaussian  quadrature  may  become  less  accurate  as  more 
points  are  used. 

C.  Behavior  of  the  Directional  Flux  in  Heterogeneous  Problems 

A  difficulty  encountered  in  transport  problems  is  that  the  flux,  as 
a  function  of  angle,  can  be  rapidly  varying  (nearly  discontinuous),  and 
that  the  discontinuities  can  exist  at  angles  that  are  not  necessarily 
parallel  to  material  interfaces.  As  an  example,  consider  the  fiat 
control  rod  in  an  otherwise  homogeneous  reactor  shown  in  figure  III-l. 
Suppose  that  for  some  energy  group,  the  neutron  mean  free  path  is  large 
compared  to  the  width  of  the  control  rod,  but  that  the  rod  is  strongly 
absorbing.  Then  the  flux  at  point  "A"  is  fairly  uniform,  except  in  those 
directions  for  which  particles  would  have  to  have  penetrated  the  rod.  In 
these  directions,  the  flux  is  depressed.  The  flux  is  rapidly  varying 
with  azimuthal  angle  for  the  direction  coming  from  the  edge  of  the 
control  rod.  Because  of  this  type  of  ill  behavior,  composite  angular 
quadrature  may  outperform  global  quadrature  in  those  kinds  of  problems 
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where  transport  methods  are  most  needed.  Some  global  quadrature  sets  are 


more  successful  in  handling  this  difficulty  than  others.  Some  of  these 
quadrature  sets  will  be  covered  next. 


Homogeneous 

Fuel-Moderator 

Region 


Fig.  III-l:  Example  Problem:  Homogeneous  Reactor  with  Control  Rod 


D.  Totally  Symmetric  Quadrature  Sets  (SN)  and  Related  Types 

Lathrop  and  Carlson  developed  a  family  of  quadrature  sets  which 

meet  various  symmetry  constraints  and  use  the  remaining  degrees  of 

freedom  to  match  angular  moments.  [Ref.  9  is  perhaps  the  most  useful  of 

many  available  references  on  the  subject.]  The  most  popular  of  these  is 

the  totally  symmetric  quadrature  sets.  These  are  the.  sets  implied  if  the 

term  "S  "  is  used  without  qualification  as  to  quadrature  set.  They 
N 

provide  symmetry  under  interchange  of  any  pair  of  axes,  under  90-degree 
rotation  about  any  of  the  three  (x,  y,  and  z)  axes,  and  under  reflection 
through  any  cardinal  (xy,  yz,  zx)  plane.  These  constraints  ensure  that 
the  results  computed  for  a  problem  will  not  depend  on  the  choice  of 
assignment  of  axis  labels,  clearly  a  physically  desirable  property. 
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However,  these  symmetries  are  of  more  importance  in  three  dimensional 


transport  than  in  the  two  dimensional  cases  of  concern  here.  The 
interchangeability  of  the  x  and  z  axes,  for  example,  is  a  waste  of 
degrees  of  freedom  that  could  better  have  been  used  matching  moments, 
for  example.  Lathrop  [Ref.  8]  has  shown  that  the  totally  symmetric 
quadratures  are  very  prone  to  ray  effect,  and  has  proposed  a 
rotationally  symmetric  quadrature  which  reduces  ray  effect.  This  set  is 
discussed  below,  along  with  other  product  quadratures. 

E.  Product  Quadrature  Sets 

Product  quadrature  sets  are  formed  as  the  product  of  independent 

A 

quadrature  sets  in  t  and  $  ,  such  that  A  =  (t  ,$  )  with  weight 

K  f  1  K  X 

* 

W,  ,  =  W  W._  and  with  k  ranging  from  1  to  K  and  1  from  1  to  L. 
Product  quadratures  of  this  type  are  particularly  appropriate  for  use  in 
two  dimensional  (xy)  geometry,  since  they  relax  the  x,z  and  y,z 
interchange  symmetries,  but  (usually)  retain  the  x,y  interchange 
symmetry,  as  discussed  above.  The  example  in  figure  III-l  showed  that 
the  flux  can  be  ill-behaved  in  azimuth.  However,  variations  in  polar 
angle  (for  any  fixed  azimuthal  direction)  result  only  in  smooth 
variations  in  flux  due  to  the  z-symmetry  of  these  problems.  In  a  sense 
then,  the  polar  quadrature  is  an  easier  problem  than  the  azimuthal 
quadrature.  Product  quadrature  sets  allow  the  apportionment  of  more 
computational  effort  to  the  latter  and  less  to  the  former.  Some  specific 
product  quadratures  are  described  below. 

*  The  set  is  usually  only  specified  for  the  principal  octant,  where  all 
three  direction  cosines  are  positive.  The  other  octants  are  defined  by 
reflections  of  the  principle  octant.  Since  only  the  x  >0  hemisphere  is 
used  in  xy-geometry,  the  total  number  of  directions  in  the  product 
quadrature  set  is  4KL. 


1.  Single  Range  Quadrature 


Single  range  quadrature  assumes  continuity  and 


differentiability  of  the  flux  throughout  the  (t,$)  plane  and  implicitly 


represents  the  flux  as  a  single  polynomial  in  t,  sin($),  and  cos($)  over 


the  entire  domain.  Lathrop  and  Carlson  developed  product  quadrature  sets 


using  Gauss-Legendre  quadrature  for  t  in  the  range  (-1,1)  and  Gauss- 


Chebyshev  quadrature  in  azimuth  [Ref.  9J.  They  reported  improved 


calculation  of  critical  radii  in  cylindrical  coordinates,  compared  to 


conventional  SN.  However,  this  scheme  is  closely  related  to  a  P^  type 


flux  representation,  and  shares  its  inability  to  exactly  represent 


vacuum  boundary  conditions.  Yvonne's  method,  DP  ,  uses  a  double  range 

N 


Gaussian  quadrature  (in  one  dimensional  problems)  to  correct  this 


deficiency.  In  two  dimensional  discrete  ordinates,  the  analogue  to  this 


approach  is  the  use  of  multiple  range  quadrature  sets. 


2.  Multiple  Range  (S„_  __)  Quadrature 

G 


Recently,  Abu-Shumays  has  developed  a  Gauss-Christof fel 


quadrature  in  sin(6)  over  the  range  (0,1)  {i.e.  0  in  (0,tt/2)}  [Ref.  2). 


This  is  a  double-range  polar  quadrature  which  takes  advantage  of  the  xy- 


plane  mirror  symmetry  of  the  problems  considered  here.  He  also  developed 


three  forms  of  quadruple-range  azimuthal  quadrature,  one  of  which 


retains  xy- interchange  symmetry.  This  azimuthal  quadrature  is  defined 


for  ipin  the  range  (0,  ir  /2)  and  extended  to  the  other  quadrants  by 


reflection  though  the  xz-  and/or  yz-plane.  It  can  exactly  represent 


vacuum  and  material  interface  boundaries  oriented  parallel  to  these 


planes.  The  combined  Gauss-Christof fel  polar  /  symmetric  quadruple-range 


azimuthal  quadrature  has  demonstrated  excellent  performance  in  reducing 


ray  effect. 
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Based  on  Abu-Shumays  findings,  discrete  ordinates  results  obtained 

with  the  combined  Gauss-Christof fel  /  symmetric  quadruple-range 

quadrature  are  here  assumed  to  be  the  best  of  the  available  Gaussian- 

type  S„  results  and  are  used  for  comparison  with  the  L  results.  The 

notation  used  to  designate  this  quadrature  is  S  indicating  K 

KG  9  LQ 

latitude  lines  with  Gauss-Christof fel  spacing  and  weights,  each  with  L 

azimuthal  points  per  quadrant  with  quadruple-range  spacing  and  weights. 

3.  Rotationally  Symmetric  Azimuthal  (S.._  .„)  Quadrature 

KG  f  LR 

Lathrop  [Ref.  8]  suggested  the  possibility  of  using  product 

quadrature  with  equally  weighted  directions  equally  spaced  in  azimuth. 

Such  a  quadrature  set  is  invariant  under  rotations  about  the  z-axis  by 

angles  which  are  multiples  of  the  azimuthal  spacing.  For  any  choice  of 

L,  this  provides  the  maximum  rotational  symmetry.  This  quadrature 

reduced  ray  effect  compared  to  conventional  sets.  Not  only  does  it 

better  model  the  rotational  symmetry  of  the  integro-dif ferential 

transport  equation,  as  Lathrop  observed,  but  it  also  has  the  advantages 

of  a  local  basis  representation  for  ill-behaved  fluxes.  This  latter 

point  is  a  consequence  of  the  fact  that  the  rotational  quadrature  set  is 

equivalent  to  a  composite  quadrature  using  the  midpoint  rule.  S.,  with 

rotational  quadrature  is  thus  a  degenerate  case  of  the  discrete  elements 

method  where  the  element  angular  quadrature  is  performed  by  the  midpoint 

rule.  The  notation  used  to  designate  discrete  ordinates  with  this 

quadrature  is  S„_  ,  indicating  K  latitude  lines  with  Gauss-Christof fel 

KG*  LR 

spacing  and  weights,  each  with  L  azimuthal  points  per  quadrant,  equally 


spaced  and  weighted. 


Lathrop  abandoned  the  rotational  quadrature  sets  with  the 
observations  that  they  reduce  but  do  not  eliminate  ray  effect  and  that 
they  require  more  directions  for  the  same  number  of  latitudes  than  the 
standard  sets:  "256  vs  144  for  N  =  16"  [Ref.  8:258].  This  implicit 
assumption  of  the  need  for  a  square  mesh  is  not  really  necessary, 
however.  It  may  be  that  only  a  few  latitudes,  with  Gauss-Christof f el 
quadrature,  will  provide  accurate  polar  quadrature,  allowing  the  use  of 
a  relatively  fine  azimuthal  mesh.  From  the  data  presented  in  reference 
8,  it  appears  that  an  S.  solution  (16  directions/octant)  has  less 

ZG  r  oR 

ray  effect  than  an  Slg  solution  (36  directions/octant)  at  44%  of  the 

cost.  An  S,  would  improve  the  accuracy  of  polar  quadrature,  with 
JGf  oR 

similar  ray  effect  and  still  have  only  67%  of  the  cost.  For  this  reason, 
SN  results  using  the  combined  Gauss-Christof fel  /  rotationally  symmetric 
quadrature  are  used  for  comparison  with  L„  results  obtained  with  the 
analogous  quadrature. 


F.  Discrete  Elements  Composite  Quadratures 

The  discrete  elements  method  uses  composite  quadratures  in  which 
the  choice  of  a  quadrature  set  entails  the  specification  of  the 
arrangement  and  weights  of  the  elements,  subject  to  the  constraint  that 
the  set  of  elements  tile  the  (t,$)  plane.  Both  one-dimensional  (slab) 
and  two  dimensional  (xy)  test  cases  are  used  in  evaluations  of  the 
discrete  elements  method.  The  quadrature  sets  used  are  described  below. 
1.  One-dimensional  Equal  Weight  (L^)  Quadrature 

In  one  space  dimension  (slab  geometry),  the  angular 
coordinates  are  oriented  so  that  the  problem  is  symmetric  in  azimuth. 
Therefore,  only  polar  quadrature  need  be  considered.  The  simplest  scheme 


is  to  divide  the  range  of  t  (-l,  +  l)  into  N  equal  intervals,  as  is 


usually  done  in  composite  numerical  integration.  The  error  bound  for  the 
quadrature  over  each  element  is  proportional  to  the  width  of  the 
interval  (i.e.  the  element  weight)  raised  to  a  power  determined  by  the 
element  quadrature  rule.  For  example,  with  Simpson's  rule  the  error  is 
proportional  to  W^.  Hence,  using  equal  weights  tends  to  minimize  the  sum 
of  the  errors  over  the  set  of  elements.  Another  advantage  is  discussed 
in  reference  to  diffusion  limit  considerations,  below. 

2.  One-dimensional  Gaussian  Weight  U-NGW>  Quadrature 

A  feature  of  the  Gauss-Legendre.  quadrature  for  one-dimensional 

S„  is  that  it.  biases  the  placement  of  directions  towards  the  poles, 

i.e.,  forward  and  backward  through  the  problem,  and  puts  reduced 

emphasis  upon  the  equatorial  (sideways)  directions.  A  composite 

quadrature  set  which  uses  elements  unequal  in  size,  such  that  the 

element  weights  are  those  of  the  Gauss-Legendre  ordinates,  would  retain 

this  feature.  It  may  be  conjectured  that  good  resolution  in  the  polar 

• 

fluxes  is  more  important  than  for  the  equatorial  fluxes,  so  that  a 
* 

"Gauss-like"  composite  quadrature  might  prove  more  effective  than  an 
equal  weight  composite  quadrature.  Gaussian  weight  composite  quadrature 
is  used  to  test  this  conjecture. 

3.  Equal  Weight  Composite  Product  (L^  ^)  Quadrature 

For  two-dimensional  (xy)  problems,  the  discrete  elements 
method,  as  derived  in  section  III-C,  uses  composite  product  quadrature. 
Although  any  collection  of  rectangular  elements  which  tiles  the  (t,$) 
coordinate  plane  could  be  employed,  a  division  into  discrete  elements  of 
equal  area  is  most  natural.  The  advantages  of  this  quadrature  are  that 
uniform  weights  minimize  the  sum  of  the  error  terms  and  that  the  set 


most  nearly  meets  the  diffusion  limit  criteria,  as  explained  below.  This 
quadrature  set  also  retains  the  advantage  of  rotational  symmetry  in 
modeling  the  rotational  symmetry  of  the  Boltzmann  transport  equation  and 
could  be  used  for  three-dimensional  problems. 


G.  Hybrid  Gauss-Christoffel/Composite  (L„„  „ )  Quadrature 

l\b  t  u 

In  two-dimehsional  problems,  as  observed  above,  the  flux  is  well- 
behaved  with  respect  to  t,  but  may  be  ill-behaved  with  respect  to  <J>. 
Consequently,  a  hybrid  quadrature  scheme  is  proposed.  This  scheme  uses 
Gauss-Christof fel  quadrature  in  the  polar  variation  and  equal  weight 
discrete  elements  for  the  azimuthal  quadrature.  The  fixed  latitudes  of 
the  polar  quadrature  take  advantage  of  the  well-behaved  character  of  the 
flux  to  provide  the  accuracy  of  a  high  order  global  quadrature,  while 
the  equal  weight  composite  quadrature  accommodates  the  potential  ill- 
behavior  of  the  flux.  This  method  provides  the  advantage  of  coupled 
angular  and  spatial  quadrature  to  steer  the  flux,  but  only  in  azimuth. 
However,  steering  in  azimuth  should  suffice  to  reduce  ray  effect.  Since 
fewer  auxiliary  fluxes  are  required,  none  being  needed  for  polar 
steering,  the  method  would  be  more  cost  effective  than  the 
quadrature,  if  it  were  at  least  as  accurate. 


K,  L 


H.  Diffusion  Limit  Considerations 


Although  numerical  transport  methods  are  needed  for  accurate  flux 
determinations  where  diffusion  theory  does  not  apply,  it  is  desirable 
that  the  numerical  method  produce  the  same  result  as  diffusion  theory 
for  problems  where  diffusion  theory  is  valid.  A  numerical  method  which 
does  this  is  said  to  satisfy  the  diffusion  limit. 
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A  discrete  ordinates  quadrature  set  which  meets  the  criteria 
enumerated  below  will  satisfy  the  diffusion  limit.  This  can  be  shown  by 
assuming  that  the  flux  is  of  the  form  ^  (£})  -  $  +  ft*  J  .  Then,  in  the 
limit  of  small  J,  discrete  quadrature  with  a  quadrature  set  meeting 
these  constraints  can  be  used  in  place  of  analytic  quadrature  in  the 
derivation  of  Pick's  law  or  of  the  diffusion  equation.  The  criteria  are: 

y  w  =  i 

L-i  m 
m 

Zy  W  =  0 
mm 
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m  m  m 
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V y  2  W  =  1/3 
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m 

yn  2  w  =1/3 
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The  physics  behind  these  constraints  is  of  interest.  Equation  (III-2) 

assures  conservation  in  that  the  calculation  of  4>  (hence  of  Q)  will  not 

lose  or  invent  particles.  Equations  (III-3)  and  (I1I-4)  serve  the  same 

purpose  with  respect  to  current.  If  the  directional  flux  is  isotropic, 

for  example,  they  ensure  that  the  current  components,  J  and  J  , 

x  y 

respectively,  are  properly  computed  as  zero.  Equation  (III-5)  causes 
cross  product  terms  to  vanish  from  the  diffusion  equation;  the  last  two 
constraints  produce  the  factor  of  1/3  in  the  diffusion  coefficient. 

In  the  limit  of  isotropic  flux,  the  flux-weighted  element  streaming 


(III-2) 

( 1 1 1— 3 ) 

(III-4) 

(III-5) 

(111-6) 
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directions  go  to  the  element  center  directions  and  the  discrete  elements 


method  becomes  equivalent  to  a  discrete  ordinates  method  with  those 

directions  and  weights  as  its  quadrature  set.  The  L  becomes,  in  this 

KG^L 

degenerate  case,  S  ,  for  example.  It  is  desirable  for  the  discrete 

KGf  LR 

element  quadrature  sets  to  meet  the  diffusion  limit. 

All  the  (xy-geometry)  quadrature  sets  described  in  this  chapter, 

and  used  in  this  report,  are  constructed  by  specifying  the  set  on  the 

primary  octant  and  filling  in  the  other  octants  by  reflection. 

Therefore,  they  have  mirror  symmetry  through  the  yz-plane,  which  ensures 

equation  (III-4)  is  satisfied,  and  through  the  xz-plane,  for  equation 

(III-3).  Either  of  these  two  symmetries,  alone,  is  sufficient  to  satisfy 

the  cross-product  constraint,  equation  (III-5).  The  conservation 

constraint  is  met  by  proper  normalization  of  the  weights.  All  of  the 

sets  used  here  are  symmetric  about  <fr  =  ir/4  and  so  are  invariant  under 

interchange  of  the  x  and  y  axes.  Thus,  if  either  of  equations  (III-6) 

and  (III-7)  are  met,  then  so  is  the  other.  The  standard  S„  quadrature 

N 

sets,  and  the  S„_  T_  sets  (for  K>1)  satisfy  enough  moment  conditions  to 
KGf  LQ 

ensure  satisfying  equations  (III-6)  and  (III  —  7) .  These  quadratures 

satisfy  the  diffusion  limit  exactly. 

An  advantage  of  the  Gauss-Christof fel  polar  quadrature  (upon  which 

Abu-Shumays  did  not  remark,  in  Ref.  2)  is  that  it  satisfies  the 

diffusion  limit  exactly,  for  K>1,  when  used  with  any  azimuthal 

quadrature  with  the  symmetries  described  in  the  previous  paragraph.  This 

2  2  2 

is  because,  in  equation  (III-6),  y  =  cos  (<f>)  sin  (0)  and  the 

m 

symmetry  about  it/4  allows  the  quadrature  points  to  be  arranged  in  pairs 

2  2 

with  complementary  azimuths.  Then  the  identity  sin  (<J>)  +  cos  (<{>)  =  1 
causes  the  choices  of  azimuths  and  weights  to  drop  out  of  equation  (III- 
6).  As  a  result,  the  Gauss-Christof fe 1  polar  /  rotationally  symmetric 
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azimuthal  quadrature,  S  and  the  corresponding  discrete  elements 

KG  t  LR 

quadrature,  L  both  satisfy  the  diffusion  limit  exactly. 

KG/L 

The  L  quadrature  (considered  as  a  fixed  quadrature  set)  does  not 
K  f  L 

meet  the  diffusion  limit  exactly.  For  the  reasons  given  in  the  last 
paragraph,  the  summation  in  equation  (III-6)  is  independent  of  the 
azimuthal  points  and  depends  only  on  the  choice  of  K.  It  reduces  to 

(1/K)  ]T  sin2  (6^  =  (1/K)  £  d-rk2)  (III-8) 

k  k 

where 

T .  =  (2k-l)  /  ( 2K)  ( III-9) 

k 

which  can  be  solved  to  obtain 

Tv  2  W  =  4  t1  +  1/ (8K2)  )  (HI-10) 

i—i  mm  J 

m 

In  practice,  equation  (III-6)  need  not  be  met  exactly,  but  the  method 

will  perform  poorly  under  some  circumstances  if  the  summation  differs 

from  1/3  significantly.  In  this  case,  the  method  converges  rapidly  to 

1/3;  for  K  as  small  as  2,  the  error  is  only  about  3%,  even  if  the 

quadrature  set  is  used  with  discrete  ordinates.  The  composite  quadrature 

of  the  discrete  elements  method  further  reduces  this  error. 

Consequently,  for  practical  purposes,  the  diffusion  limit  is  met. 

Similar  arguments  can  be  made  for  the  one-dimensional  quadrature,  L^, 

1  2 

for  which  the  summation  is  given  by  —  (1  -  1/N  )  Also,  it  can  be 
shown  that  these  equal  weight  schemes  are  closest  to  meeting  the 
diffusion  limit,  in  that,  if  the  weight  of  one  element  is  increased  and 
that  of  another  correspondingly  decreased,  the  approximation  to  1/3  is 
less  accurate. 
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I .  Summary 


For  xy-geometry,  there  are  two  natural  choices  of  quadrature  for 

the  discrete  elements  method.  With  steering  of  the  streaming  angles  in 

both  polar  and  azimuthal  coordinates,  the  equal  weight  product 

quadrature,  L  is  used.  With  steering  only  in  azimuth,  the  hybrid 
K,  L 

Gauss-Christof fel  polar  /  equal  weight  (rotational ly  symmetric) 

azimuthal  product  quadrature,  L  is  used. 

KG  /  L 

The  performance  of  discrete  ordinates  is  strongly  dependent  upon 
the  choice  of  quadrature  set,  and  the  optimum  choice  is  problem 
dependent.  A  meaningful  performance  comparison  between  discrete  elements 
and  discrete  ordinates  requires  comparison  with  several  discrete 
elements  quadratures: 

1  -  Conventional  (Totally  Symmetric)  SN 

This  is  often  a  sub-optimal  performer,  but  the  production 

codes  in  current  use  are  designed  around  this  quadrature  set,  so  that  it 

is  used  by  default.  Product  quadratures  are  awkward  to  use  and 

inefficiently  implemented.  DOT-4.3  is  an  example  [Ref.  12].  In  a 

practical  sense,  S„  is  the  standard  of  performance. 

N 

2  -  Gauss-Christof fel  /  Rotationally  Symmetric  S 

KG  f  LR 

This  quadrature  underlies  the  hybrid  L  scheme. 

KG  t  Li 

Comparison  of  these  two  will  determine,  if  the  discrete  elements  method 
works  well,  whether  it  works  because  of  its  underlying  quadrature  set  or 
because  of  its  composite  quadrature  within  each  element. 


3  -  Equal  Weight  /  Rotationally  Symmetric  S 

KE  r  LR 

This  quadrature  underlies  the  L  scheme  and  comparison 

K/L 

would  have  the  same  benefits  as  in  the  previous  paragraph.  However, 
SKE,LR  never  as  accurate  as  S^,  LR  while  showing  essentially  the  same 
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ray  effects.  These  two  Sv,  methods  have  equal  computational  cost. 

N 

Therefore,  comparison  of  L„  with  S„_  is  reported  here. 

L  KG»LR 

4  -  Gauss-Christof fel  /  Quadruple  Range  S„_ 

KG  t  LQ 

This  is  a  special  quadrature  set,  ideally  adapted  to  the 
test  problems  used  here.  Comparison  with  this  quadrature  is  a  demanding 
test  of  the  discrete  elements  method.  It  should  be  noted,  however,  that 
such  ideal  discrete  ordinates  quadratures  are  not  available  for  many 
practical  problems.  For  example,  the  presence  of  diagonally  oriented 
absorbers,  material  interfaces,  or  vacuum  ducts  would  invalidate  the 
assumptions  of  the  quadruple-range  azimuthal  quadrature  but  would  not  be 
expected  to  degrade  the  performance  of  discrete  elements. 

All  the  quadrature  schemes  considered  here  meet,  or  very  nearly 
meet,  the  diffusion  limit  constraints,  equations  (III-2)  through 


IV.  Quadrature  Rules  for  Discrete  Element  Flux  Model line 


The  previous  chapter  noted  that  there  are  many  possible  choices  of 

quadrature  sets  for  discrete  ordinates.  The  corresponding  freedom  in  the 

discrete  elements  method  is  in  the  choice  of  quadrature  rule  within  each 

element.  This  quadrature  rule  determines  several  characteristics  of  the 

flux  model:  whether  it  is  piecewise  continuous  or  discontinuous,  the 

order  of  polynomial  used  as  the  basis  function  within  each  element,  and 

the  possible  range  of  variation  of  the  flux-weighted  streaming 

direction.  Three  quadrature  rules  were  tested:  Simpson's  rule,  three- 

point  Gauss-Legendre  rule,  and  the  Newton-Cotes  (five-point)  rule.  This 

chapter  reviews  the  rules  and  their  Taylor  series  error  bounds  and 

considers  their  relative  advantages  and  disadvantages.  Table  IV-1 

provides  the  formulas  for  flux-weighted  mean  angle,  optimized  for 

minimum  computational  operations.  The  parameter  "s"  in  these  formulas 

represents  T  or  <J> .  The  mean  azimuth  formulas  are  for  use  with  L  and 

K,L 

L  T.  Corresponding  mean  t  formulas  are  used  with  L  and  L  . 

K0t,JL.  K,L  N 


A.  Simpson's  Rule 

Simpson's  (three-point)  rule  is  the  lowest-order  scheme  considered. 
A  convenient  form  for  use  here  is 


/  f (x)  dx  =  h  (f(a-h)  +  4  f(a)  +  f(a+h)]  /  3 
a-h 


(IV-1) 


An  error  bound  (Taylor  series  residual,  R)  for  this  approximation  is 


5  (4) 

R  =-  0.000347  (Ax)  fv  '(x’) 


(IV-2) 


where  Ax  =  2h  and  x'  is  some  point  in  (a-h,a+h). 


v  v  v  v  >  ■. 
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With  Simpson’s  rule,  the  directional  flux  is  modelled  by  a 
continuous  piecewise-quadratic  polynomial.  Continuity  is  accomplished 
not  by  overlapping  finite  elements,  but  implicitly,  through  the  use  of 
the  same  auxiliary  flux  at  the  boundary  of  adjacent  discrete  elements. 
If  the  actual  flux  is  discontinuous,  or  nearly  so,  then  the  model  will 
be  poor  only  in  the  element(s)  in  which  the  discontinuity (s)  occur. 
Discontinuities  at  the  boundaries  between  (directional)  octants  are 
included  in  this  model,  since  the  auxiliary  fluxes  along  octant 
boundaries  are  independently  computed  for  each  octant  as  part  of  the 
spatial  mesh  walk.  Consequently,  vacuum  boundaries  parallel  to  principal 
planes  can  be  modelled  exactly  by  the  discrete  elements  method, 
regardless  of  the  choice  of  quadrature  rule. 

Advantages  of  Simpson's  rule  include: 

-  Mean  streaming  direction  can  vary  over  full  range  of  the 
element.  This  could  be  valuable  in  problems  with  streaming  ducts  through 
shields  since  particles  could  be  modelled  as  streaming  directly  down  the 
duct,  regardless  of  its  orientation. 

-  Some  auxiliary  directions  are  on  boundary  between  elements. 
These  can  be  used  twice,  if  sufficient  storage  is  available,  reducing 
computational  cost  of  auxiliary  fluxes. 

B.  Three-Point  Gauss-Legendre  Rule 

The  Gauss-Legendre  (three-point)  rule  is  of  higher  order  than 
Simpson's  rule,  but  has  the  same  computational  cost  (assuming  auxiliary 
fluxes  are  not  reused  due  to  storage  constraints).  It  is  the  only  method 
used  here  which  models  the  flux  as  discontinuous  piecewise  polynomial.  A 
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convenient  form  for  use  here  is 


a+h 

f  f(x)  dx  =  h  [5  f(a-bh)  +  8  f(a)  +  5  f(a+bh)]  /  9  (IV-3) 

a-h 

1/2 

where  the  auxiliary  flux  offset  factor  is  b  =  (3/5) 

An  error  bound  (Taylor  series  residual,  R)  for  this  approximation 
is 

R  =  -  4.96xl0-7  (Ax)7  f(6)(x')  (IV-4) 

where  Ax  =  2h  and  x'  is  some  point  in  (a-h, a+h) . 

Advantages  of  the  three-point  Gauss-Legendre  quadrature  are: 

-  Higher  order  on  individual  discrete  elements  with  only 
slightly  higher  cost  (no  more  cost  than  Simpson's  rule  if  boundary 
auxiliary  fluxes  aren't  reused  due  to  memory  constraints) 

-  Broad  range  of  variation  of  mean  angle,  compared  to  two- 
point  Gauss-Legendre,  but  mean  angle  cannot  reach  the  edge  of  the 
element.  This  could  interact  well  with  spatial  quadrature  schemes  since 
they  are  inaccurate  for  very  small  y  or  h . 

-  The  underlying  model  is  equivalent  to  a  piecewise- 
discontinuous,  piecewise-quartic  (fourth-order)  polynomial  fit  to  the 
directional  flux.  In  this  model,  if  any  discontinuities  occur  at  (or 
very  near)  to  the  boundary  between  elements,  they  may  be  modelled 
accurately,  which  is  a  potential  advantage.  This  provides  a  simple  way 
to  optimize  the  method  for  problems  where  such  discontinuities  are 
anticipated  to  occur  along  a  limited  number  of  known  directions. 
(Examples:  hexagonally  shaped  or  diagonally  oriented  absorbers,  sources, 
or  material  interfaces  in  an  otherwise  rectangular  system.) 

IV-3 


■c. 
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C.  Newton-Cotes  Closed  (Five-Point)  Rule 

The  Newton-Cotes  (five-point)  rule  is  of  the  same  order  as  the 
Gauss-Legendre  (three-point)  rule,  and  the  coefficients  of  the  error 
bounds  are  approximately  equal.  Newton-Cotes  rule  provides  equivalent 
accuracy,  but  at  greater  computational  cost  due  to  the  increased  number 
of  auxiliary  fluxes  needed.  It  is  tested  for  use  in  the  discrete 
elements  method  because  it  provides  continuity  at  the  boundaries  between 
elements,  unlike  the  Gauss-Legendre  rule.  A  convenient  form  for  use  here 
is 

a+h 

/  f (x)  dx  =  h  (7  f(a-h)  +  32  f(a-h/2)  +  12  f(a) 
a-h 

+  32  f (a+h/2)  +  7  f(a+h)]  /  45  (IV-5) 

An  error  bound  (Taylor  series  residual,  R)  for  this  approximation  is 

R  -  -  5.17xl0~7  (Ax)7  f(6)(x')  ( IV-6) 

where  Ax  =  2h  and  x'  is  some  point  in  (a-h, a+h). 

The  advantages  are: 

-  Computational  cost  of  using  five  auxiliary  directions  is 
offset  by  the  fact  that  two  of  them  are  on  the  element  boundaries  and 
can  be  reused,  as  with  the  Simpson's  rule. 

-  Mean  streaming  direction  can  vary  over  the  full  range  of  the 
discrete  element,  as  with  Simpson's  rule,  which  could  be  beneficial  in 
duct  streaming  problems,  for  example 

-  The  underlying  model  is  a  continuous  piecewise-quartic  fit 
to  the  directional  flux.  The  comments  about  the  Simpson's  rule  model 


apply  here,  as  well. 


Simpson's  Rule: 


<s>  =  sQ  +  h 


f(sQ+h)  -  f(sQ-h) 


f(s  +h)  +  4  f(s  )  +  f(s  -h) 


( IV- 7) 


where  h  =  As/2 


Gauss-Legendre  Three-Point  Rule: 


<s>  =  s0  +  g 


f(sQ+g)  -  f (sQ-g) 


f(sQ+g)  +1.6  f(sQ)  +  f(sQ-g) 


( IV-8) 


where  h  =  As/2 


g  =  (3/5) 1/2  h 


Newton-Cotes  Five-Point  Rule: 


<s>  =  sQ  +  h 


f(sQ+h)  -  f(s0-h)  +  c1  [f (sQ+h/2)  -  f (sQ-h/2) ] 
f(s0+h)  +  f(sQ-h)  +  c2  I f (Sg+h/2)  +  f (sQ-h/2) ]  +  c3  f(sQ) 


(IV-9) 


where  h  =  As/2 
cx  =  16/7 
c2  =  32/7 
c3  -  12/7 


Table  IV-1:  Flux-Weighted  Mean  Angle  Formulas 


IV-5 


D.  Summar 


The  discrete  element  quadrature  rule  is  used  in  computing  the  flux- 
weighted  mean  streaming  direction  for  the  element.  Its  only  effect  on 
the  method  is  through  this  steering  of  particle  streaming.  Three 
quadrature  rules  were  selected  for  evaluation.  Simpson’s  rule  is  low 
order  but  computationally  inexpensive  and  allows  the  streaming  direction 
to  be  steered  anywhere  within  the  element  of  angle.  With  the  same 
expense  (three  auxiliary  flux  calculations  per  element),  Gauss-Legendre 
provides  higher  order  and  the  possible  advantages  of  a  piecewise 
discontinuous,  piecewise  polynomial  flux  representation.  Gauss-Legendre 
restricts  the  steering  of  the  streaming  angle  to  the  center  77%  of  the 
element  range  (in  one  angular  dimension),  which  may  limit  its  ability  to 
model  a  streaming  duct,  for  example,  but  also  avoids  space  quadrature  in 
directions  nearly  parallel  to  the  sides  of  the  space  cell,  for  which 
most  space  quadrature  schemes  have  very  poor  accuracy.  The  Newton-Cotes 
rule  uses  five  auxiliary  fluxes  per  element  and  so  is  more  expensive  to 
compute.  It  uses  a  continuous,  piecewise  polynomial  flux  representation 
and  allows  the  streaming  direction  to  be  steered  anywhere  in  the 
element. 

The  above  discussions  of  error  bounds  and  accuracy  all  implicitly 
assume  that  the  auxiliary  fluxes  (f(sQ),  etc.)  are  exact.  However,  they 
are  calculated  by  numerical  means,  discrete  ordinates  spatial 
quadrature.  The  errors  in  the  data  for  the  quadrature  rules,  as  it  were, 
may  well  dominate  their  performance.  The  results  of  numerical 
experimentation  are  reported  in  later  chapters.  The  spatial  quadrature 
methods  employed  here  are  discussed  in  the  next  chapter. 


V.  Spatial  Quadrature  Methods 

The  research  reported  in  this  dissertation  is  primarily  concerned 
with  two-dimensional  Cartesian  (xy)  geometry,  although  the  special  case 
of  one-dimensional  Cartesian  (slab)  geometry  is  also  considered.  The 
discrete  element  method  uses  the  same  spatial  integration  methods  as  the 
discrete  ordinates  method;  therefore,  these  spatial  methods  are  reviewed 
in  this  chapter. 

A.  One-Dimensional  (Slab)  Spatial  Quadrature 

Four  methods  of  one-dimensional  spatial  quadrature  were  evaluated 
in  this  study;  step,  diamond  difference  (with  negative  flux  fixup),  step 
characteristic,  and  linear  characteristic.  Several  recent  papers  have 
concerned  these  methods.  Alcouffe,  et  al.  [Ref.  3],  proposed  the  linear 
characteristic  (L.C.)  method  and  reported  excellent  performance,  in 
comparison  with  other  methods,  based  on  numerical  testing.  Larsen  and 
Miller  [Ref.  6]  evaluated  the  convergence  (vs.  space  cell  size)  of  these 
methods,  as  did  Lee  and  Vaidyanathan  [Ref.  11).  The  concepts  behind  the 
methods,  their  formulas,  and  some  of  their  advantages  and  disadvantages 
are  summarized  below. 

1.  Coordinates  and  Symbols 

In  the  one-dimensional  coordinate  system,  x  is  the  space 
coordinate;  the  angular  coordinates  are  rotated  (compared  to  xy 


geometry)  such  that  0  is  the  angle  with  respect  to  x,  and  H  is  the  x 


direction  cosine  (  y  =  cos{8)).  For  convenience  in  computer 


implementation,  let  be  the  directional  flux  entering  the  space  unit 

cell  of  width  Ax  and  F  be  that  leaving.  This  convention  allows  the 

out 

use  of  the  absolute  value  of  y  and  requires  that  the  program  explicitly 

keep  track  of  whether  F.  represents  F.  _  and  F  .  represents  F  .  .. 
c  in  left  out  right 

(for  positive  y,  i.e.,  for  right-bound  particles)  or  vice  versa  (for 
negative  y,  i.e.,  for  left-bound  particles).  Using  this  scheme,  the 
spatial  quadrature  method  need  only  consider  positive  directions. 
Negative  directions  are  handled  as  the  reflection  of  positive  ones.  With 
this  notation,  the  transport  equation  to  be  solved  is: 


y  dF/dx  +  0  F  =  og 

(V-l) 

0-0  total 

(V-2) 

c  =  a  ..  /  a  ^  , 

scatter  total 

(V-3) 

Q  =  c  $  +  S/a 

(V-4) 

The  effective  source  parameter,  Q,  has  dimensions  of  a  directional  flux, 
and  is  physically  significant  as  the  value  the  flux  asymptotically 
approaches  as  y  0.  Defining  Q  in  this  way  is  computationally 
convenient  and  efficient,  as  well  as  physically  meaningful,  and  is 
consistent  with  the  usage  in  chapter  II. 


Balance  Equation 


The  balance  equation  is  a  conservation  relation  among  the 


cell-edge  and  cell-average  fluxes  obtained  by  integrating  equation  (V-l) 


across  a  space  cell: 


(V-5) 


y  (F 


out 


-  Fin> 


/A  x 


0  F 


o  Q, 


where  Fc  is  the  cell  average  flux  and  o  Qc  is  the  cell  average  effective 
source.  (The  subscript  on  Q  will  be  omitted,  except  where  necessary  to 
avoid  ambiguity.)  The  cell  optical  thickness,  t  ,  is  the  parameter 
defined  as 


t  =  aA  x  /y 


( V— 6 ) 


Using  this  parameter,  the  balance  equation  becomes 


(F 


out 


F. 

in 


)  /e 


+  F  =  Q 
c 


( V— 7 ) 


3.  Conservation  Considerations 

Conservation  of  particles  is  assured  if  two  conditions  are 

met.  First,  the  spatial  quadrature  scheme  must  be  based  upon  the  (exact) 

balance  equation,  regardless  of  any  approximate  auxiliary  relations  it 

assumes.  This  condition  ensures  conservation  of  particles  inside  each 

spatial  mesh  cell.  All  the  spatial  quadratures  used  here  meet  this 

condition.  Secondly,  particles  must  be  conserved  in  crossing  the 

interfaces  between  cells.  In  the  discrete  ordinates  method,  this  is 

achieved  by  using  the  flux  out  of  one  cell  as  the  flux  into  the  next 

cell.  More  precisely,  it  is  the  component  of  current  in  the  direction 

normal  to  the  face  of  the  cell  which  must  be  carried  unchanged  across 

the  cell  boundaries.  In  one  dimension,  this  current  is  J  =  y  F  .  In 

x  mm 

the  discrete  ordinates  method,  the  directions,  y  ,  are  fixed  so  that 

m 

continuity  of  flux  across  the  cell  interfaces  is  sufficient. 
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In  the  discrete  elements  method,  this  same  logic  applies  to  the 
auxiliary  fluxes,  whose  directions  are  fixed.  However,  the  main  fluxes 


in  the  discrete  element  method  use  steered  streaming  directions  which 
are  not  fixed  (which  is  the  major  objective  of  discrete  elements). 
Conservation  is  accomplished  in  the  same  two  steps.  Within  each  space 


cell,  the  streaming  angle  for  each  element  is  assumed  to  be  constant.  It 
may  have  a  different  value  after  each  iteration,  but  while  the  spatial 
quadrature  is  performed,  the  flux-weighted  streaming  directions  are 


treated  as  fixed.  This  assures  conservation  within  each  space  cell, 
since  the  explicitly  conservative  spatial  quadrature  methods  are  used, 
as  in  S  .  Conservation  across  cell  interfaces  is  achieved  by  explicitly 
using  the  normal  current  out  of  one  cell  as  the  normal  current  into  the 
next  cell.  Since  the  streaming  directions  are  discontinuous  across  cell 
interfaces,  so  are  the  fluxes,  but  in  just  such  fashion  as  to  maintain 


continuity  of  the  current.  Consequently,  there  are  two  physically 
meaningful  results  computed  by  the  method,  the  cell-average  fluxes  and 
the  cell-boundary  currents. 


4.  Step  Method 

The  step  method  assumes  that  the  flux  is  constant  across  each 
space  cell,  with  discontinuities  at  the  inbound  edges.  This  is 
equivalent  to  assuming  the  auxiliary  relation: 


F  =  F 
,c  out 


(V-8) 


This  leads  to  a  solution  of  the  balance  equation: 


Fc-Fout=2  e/  (1+e)  +  Pin'  (1+e) 


(V-9) 
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This  solution  has  several  computational  advantages.  The 
coefficients  of  Q  and  are  in  the  range  (0,1)  for  all  possible 
optical  thicknesses.  Therefore,  the  numerical  method  which  uses  these 
coefficients  is  stable  and  absolutely  convergent.  The  method  is 
positive,  in  that,  for  any  positive  inputs  (Q  and  F^  ) ,  the  computed 
fluxes  are  positive;  hence,  no  negative  flux  "fixup"  is  required. 
Positive  methods  also  have  the  benefit  of  producing  solutions  which  vary 
smoothly  from  cell  to  cell  across  the  space  mesh.  Finally,  there  are  no 
special  cases,  so  that  no  if-tests  are  required.  This  combines  with  the 
computational  simplicity  of  equation  (V-9)  to  make  the  method  extremely 
fast. 

The  disadvantage  of  this  method  is  its  poor  accuracy,  which 
requires  a  prohibitively  fine  space  mesh  to  provide  usable  results  for 
most  purposes.  Since  the  method  is  smooth  and  positive,  it  may  be 
adequate  for  computing  the  auxiliary  fluxes  used  in  steering  the 
streaming  directions.  For  that  application,  only  relative  levels  are 
important,  as  is  seen  from  the  form  of  equations  (11-20)  and  (11-21). 

5 .  Diamond  Difference 

The  diamond  difference  method  assumes  that  the  flux  varies 
linearly  across  the  space  cell.  Thus,  the  cell-average  flux  is 
approximated  by  the  average  of  the  cell-edge  fluxes.  This  is  equivalent 
to  the  familiar  "diamond  difference"  approximation  for  numerical 
differentiation  and  gives  the  auxiliary  relation: 


(F  +  F.  )  /  2 
out  in 


(V-10) 


-V\-Vv 


This  combines  with  the  balance  equation  to  give  a  solution: 


F  =  Q  2  e / (2  +  e  )  +  F.  (2-e)/(2+e)  (V-ll) 
out  m 

F^  is  then  computed  directly  from  the  auxiliary  equation. 

This  method  has  several  advantages.  It  is  numerically  stable  and 
absolutely  convergent;  the  coefficients  of  F  and  Q  in  the  explicit 
solution  for  F^  lie  in  the  range  (-1,+1),  since  the  optical  thickness  is 
positive.  As  in  the  step  method,  the  calculations  are  simple  and  no  if- 
tests  are  required.  This  method  is  substantially  more  accurate  than  the 
step  method,  and  is  easily  generalized  to  curvilinear  geometries,  as 
with  cylindrical  or  spherical  coordinates. 

The  diamond  difference  method,  however,  is  not  a  positive  method. 
It  may  predict  negative  fluxes,  which  are  not  physically  reasonable. 
Relatively  fine  space  meshes  are  needed  to  avoid  negative  fluxes.  The 
calculation  of  space-integral  quantities,  such  as  total  absorption 
rates,  is  insensitive  to  the  presence  of  a  few  negative  fluxes;  but  the 
pointwise  values,  such  as  scalar  flux  as  a  function  of  position,  tend  to 
be  oscillatory  and  inaccurate.  These  inaccuracies  impair  its  usefulness 
for  computing  the  auxiliary  fluxes,  but  it  may  be  adequate  for  the  main 
fluxes  of  the  discrete  elements  method. 

6.  Negative  Flux  Fixup  for  the  Diamond  Difference  Method 

In  order  to  avoid  the  (conceptually  disconcerting)  negative 
fluxes,  a  fixup  scheme  is  often  employed.  This  is  done  by  computing  FQUt 
using  the  normal  diamond  difference  equations,  as  shown  above.  If  the 
result  is  negative,  however,  it  is  (arbitrarily)  set  to  zero.  This  new 
auxiliary  relation  (F  .  =  0)  is  substituted  into  the  balance  equation, 


which  is  then  solved  for  F  : 

c 

Fc  =  Q  +  F.n  /e  (V-12) 

The  advantage  of  this  scheme  is  that  the  pointwise  fluxes  are  more 
reasonable,  while  the  space-integral  values  are  nearly  as  accurate  as 
before,  assuming  the  fixup  is  invoked  only  infrequently. 

The  disadvantage  of  the  fixup  is  that  the  numerical  method  becomes 
nonlinear  and  convergence  is  no  longer  assured.  An  oscillatory 
instability  can  result  wherein  the  fixup  is  required  on  one  iteration, 
but  the  fixed-up  flux  avoids  the  need  for  the  fixup  on  the  next 
iteration,  yet  that  un-fixed-up  flux  requires  the  fixup  on  the 
subsequent  iteration,  etc.  This  can  happen  if  the  fixup  is  required  too 
frequently,  which  is  sometimes  the  case  for  realistic  problems. 

7.  Considerations  Regarding  the  Use  of  Fixups 

There  are  two  schools  of  thought  regarding  the  use  of  negative 
flux  fixups.  One  viewpoint  is  that  the  fixup  degrades  integral  values 
and  should  not  be  used.  Negative  fluxes  are  seen  as  being  merely 
inaccurate  numerical  representations  of  the  (correct)  positive  values. 
Further,  the  presence  of  too  many  negative  fluxes  is  considered  an 
indication  of  a  need  to  refine  the  space  mesh.  The  other  viewpoint  is 
that  negative  flux  fixups  should  be  used  since  pointwise  values  are 
improved,  so  long  as  integral  values  are  not  significantly  altered  by 
the  fixup.  Changes  in  integral  values  are  seen  as  an  indication  of  the 
need  to  refine  the  space  mesh.  Since  amelioration  of  ray  effects  is  an 
objective  of  this  research,  pointwise  values  are  of  importance,  and  the 
latter  point  of  view  is  taken  here. 


8.  Step  Characteristic  Method 


The  step  characteristic  method  is  derived  by  assuming 
streaming  along  the  "characteristic  lines"  of  the  transport  differential 
equation,  i.e.,  in  direction  .  The  transport  equation  (V-l)  is 
integrated  analytically,  with  the  assumption  that  Q(x)  is  a  constant, 
Q^,  across  the  space  cell.  This  provides  an  auxiliary  equation  for  F  : 


F 

out 


e"E  F.  +  (1  -  e_E)  Q 

in  * 


(V-13) 


The  resulting  value  is  then  used  in  the  balance  equation  to  obtain 


F 

c 


Q  + 


(F.  - 

in 


F  J 

out 


/  e 


(V-14) 


An  advantage  of  this  method  is  its  positivity,  which  avoids  the 
need  for  negative  flux  fixups.  The  method  also  avoids  the  point-to-point 
oscillations  of  the  diamond  difference  method.  With  both  smoothness  and 
accuracy,  the  method  should  be  applicable  to  both  the  auxiliary  and  main 
fluxes  of  the  discrete  elements  method. 

The  disadvantage  of  the  method  is  the  computational  cost  of 

evaluating  the  exponential  function.  This  cost  can  be  minimized  if  the 

s 

coefficients  are  computed  only  once  and  stored  in  arrays,  assuming 
sufficient  storage  is  available  for  this  use,  -  or  are  computed  only  upon 
changing  material  regions,  being  passed  to  the  step  characteristic 
subroutine  as  a  parameter.  * 

It  is  interesting  that  (in  this  one-dimensional  case,  only)  the 
diamond  difference  method  is  equivalent  to  the  step  characteristic 
method  with  the  exponential  approximated  by  (2-e)/(2+e).  The  leading 
error  term  in  this  approximation  is  e  ^/12. 


9.  Linear  Characteristic  Method 


i 

I 


The  linear  characteristic  method  is  derived  as  was  the  step 

characteristic  method  but  with  the  improved  approximation  that  the 

effective  source,  Q(x),  is  assumed  to  vary  linearly  across  the  space 

cell.  In  order  to  estimate  this  variation,  it  is  necessary  to 

accumulate,  in  the  outer  iteration,  not  only  the  cell  average  effective 

source,  Q  ,  but  also  the  edge  values  for  each  cell,  Q.  and  Q  .  (Note 
*c  ^  *in  **out 

that  O  of  one  cell  and  Q.  of  the  next  cell  are  different  if  the 
*out  m 

cross-sections  are  different.)  Then  Q(x)  can  be  approximated  by 


Q(x)  =  Qc  +  P  (x-xc)  /Ax 


where 


P  =  Q. 


out 


in 


(V-15) 


(V-16) 


However,  if  jpj  >  2  Qc  ,  then  Q{x)  will  be  negative  at  one  end  of  the 
cell.  Therefore,  a  negative  source  fixup  is  used,  which  reduces  the 
slope  of  Q(x)  so  that  it  just  reaches  zero  at  the  edge  of  the  cell.  This 
is  done  by  recomputing  P  (if  |p|  >  2  Qc)  as 

P  =  2  Q  SIGN (Q  -Q.  )  (V-17) 

out  in 

The  transport  equation  is  integrated  analytically  to  obtain  FQUt  in 
terms  of  F.  ,  Q,  P,  and  the  optical  thickness,  e  : 


% 

V 

V 


V 


F  =  F.  e_e+  Q  (1  -  e"e)  +  P  [1  -  (1/2  +  1/e)  (1  -  e"E)  ]  (V-18) 

out  in 


This  result  is  used  in  the  balance  equation  to  obtain  F 


=  Q  +  (Fin  "  Fout)  /l 


V-9 
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Alcouffe,  et  al.,  have  reported  that  the  advantage  to  this  method 
is  that  it  is  very  much  more  accurate  than  the  other  methods  considered. 
Although  it  requires  more  storage  and  computation  per  cell,  it  provides 
equivalent  accuracy  with  a  much  coarser  mesh,  so  that  storage  and 
execution  costs  are  both  reduced  for  a  given  accuracy  of  results 
[Ref.  3].  The  disadvantage  of  somewhat  more  complicated  programming  is  a 
minor  one.  This  method  is  appropriate  to  both  the  auxiliary  and  main 
fluxes  of  the  discrete  element  method. 


B.  Two-Dimensional  (xy)  Spatial  Quadrature 


Three  methods  of  xy  spatial  quadrature  are  considered:  step, 
diamond  difference  (with  negative  flux  fixup),  and  step  characteristic. 
Fine  mesh  spacing  is  used  with  these  methods  to  obtain  adequate 
accuracy,  rather  than  the  higher  order  methods  which  have  been  presented 
very  recently.  Conservation  is  achieved  by  the  same  means  as  in  one 
dimension,  as  described  above. 

1 .  Coordinates  and  Symbols 

The  two-dimensional  space  and  angular  coordinate  systems  used 
to  define  the  angular  quadratures  are  used  here  for  the  spatial 
quadratures.  The  unit  space  cell  is  shown  in  figure  V-l.  As  in  the  one¬ 
dimensional  case,  the  direction  cosines  are  all  assumed  positive  and  the 
other  three  octants  are  obtained  by  reflections.  Due  to  the  z-symmetry, 
only  the  upper  hemisphere  of  directions  (  x  >  0)  need  be  considered.  The 


fluxes  entering  the  cell  are  F.  and  F.  and  are  the  fluxes  on  the 

in,x  in,y 

left  ^x=xi_j/2^  and  bottom  ^y=yi-i/2*  faces»  respectively.  Similarly, 

F  and  F  .  are  the  fluxes  on  the  right  and  top  faces  of  the  cell. 

out,x  out,y 

A 

These  fluxes  are  directional  fluxes  in  a  fixed  quadrature  direction,  n. 


2.  Balance  Equation 


Using  the  above  notation,  the  transport  equation  becomes 


UdF/dx  +  t\  dF/dy  +  a  F  =  a  Q  (V-20) 

where  Q  is  defined  by  equation  (V-4),  as  before.  Integrating  over  the 
area  of  the  space  cell  yields  the  balance  equation: 


(F  -  F.  )  /a 

out,x  m,x 


(F 


out,y 


-  F.  )  /B 
m,y 


=  2, 


(V-21) 


where  the  optical  thicknesses  are 


a  =  a  A  x  /v  (V-22) 

8  =  a  A  y  /n  (V-23) 


and  the  cell-edge  fluxes  are  averages  over  their  respective  cell  faces, 

F  is  the  cell-average  flux,  and  Q  is  the  cell-average  source  term.  The 
c  c 

subscript  on  Q  will  be  dropped  except  where  required  to  avoid  ambiguity. 

3 .  Step  Method 

As  in  slab  geometry,  the  step  method  in  xy  geometry  assumes 
auxiliary  relations 


p  =  F  =  p 

out,x  out,y  c 


(V-24) 


This  assumption  reduces  the  balance  equation  to  a  single  unknown.  The 


balance  equation  is  solved  to  obtain 


F 

c 


Q  +  (1/a)  F.  +  (1/6;  F. 

_ _ m,x _ m,y 


1  +  (l/o) 


(1/  B) 


(V-25) 


The  step  method  has  the  same  advantages  and  disadvantages  in  two 
dimensions  as  were  described  above  for  the  one-dimensional  case.  Another 
advantage  is  that  it  can  be  readily  extended  to  curvilinear  coordinates. 

4.  Diamond  Difference 

The  diamond  difference  method  assumes  that  the  flux  is  linear 
in  both  x  and  y  over  a  space  cell,  so  that  the  cell-average  flux  is 
equal  to  the  arithmetic  average  of  the  cell-edge  fluxes  on  opposite 
sides  of  the  cell.  As  before,  the  cell-edge  fluxes  represent  averages 
over  their  respective  faces  of  the  cell.  This  assumption  leads  to  an 
auxiliary  relation  for  each  dimension: 


F 

c 

=  (F.  +  F  )  /  2 

m,x  out,x 

(V-26) 

F 

c 

=  (F.  +  F  )  /  2 

m,y  out,y 

(V-27) 

These  relations,  together 

with  the  balance  equation  form 

an  algebraic 

system  which  is  solved  for 

F  t  . 
out,x 

F  ,  =  2  ♦ 

F,  (a  -  b  -  1/2)  +  F.  2b 

in,x  m,y 

(V-28) 

out  f  X 

a  +  b  +  1/2 

where,  for  computational  efficiency,  a  =  1/a  and  b  =  1/6  . 

The  remaining  fluxes,  F  and  F  ^  ,  are  then  evaluated  using  the 

c  out,y 


auxiliary  relations. 


The  diamond  difference  method  provides  good  accuracy  for  global, 
i.e.  space-integrated,  quantities,  such  as  total  absorption;  but  the 
method  is  not  positive  and  predicts  negative  fluxes  if  the  space  mesh  is 


too  coarse.  Lathrop  has  researched  the  issue  of  positivity  versus 
accuracy.  He  states  that: 

"A  version  of  the  variable  weight  scheme  in  which  weights 
depend  not  only  on  the  space-angle  mesh  but  also  on  particle 
sources  and  fluxes  is  suggested  as  a  means  of  obtaining  the 
highest  accuracy  consistent  with  a  positive  difference  scheme, 
but  it  is  noted  that  such  schemes  are  computationally  more 
expensive  than  available  corrective  recipes  used  in 
conjunction  with  nonpositive  schemes."  [Ref.  7:475] 

The  only  variable  weight  scheme  used  here  is  the  step  characteristic 

method.  It  will  be  discussed  after  the  corrective  recipes,  which  are 

considered  next. 


5.  Negative  Flux  Fixup  for  the  Diamond  Difference  Method 

The  design  and  application  of  negative  flux  fixup  in  two 
dimensions  is  similar  to  that  used  in  one  dimension,  with  the  added 
complications  that  one  or  the  other  or  both  of  the  cell-exiting  fluxes 
can  be  negative,  so  several  special  cases  must  be  handled.  The 
algorithm,  optimized  for  minimum  if-tests  and  computer  operations, 
proceeds  as  follows: 


1  -  Compute  F  ^  from  equation  (V-28) 

out,x 

2  -  If  F  ^  is  not  negative,  compute  F  ^  from  the  diamond 

out,x  out,y 

difference  relations: 


a  -  If  F  is  not  negative,  compute  F  from  either  diamond 

out,y  c 

difference  equation,  which  completes  the  calculations. 

b  -  If  F  is  negative  (but  F  ^  was  not),  then  replace 

out,y  out,x 

the  y  diamond  relation,  equation  (V-27),  with  the  fixup  Fout  y  =  0 


and  solve  the  resulting  system  of  equations  for  F 


out,x 


„  Q  +  F.  (a  -  1/2)  +  F.  b 

F  =  in,x  in,y 

out,x  — — - - - 

a  +  1/2 


(V-30) 


If  the  revised  F  is  not  negative,  compute  F  from  the  x  diamond 

out,x  c 

relation,  equation  (V-26),  which  completes  the  calculations.  If  the 

revised  E'out  x  is  negative,  then  both  outgoing  fluxes  require  fixups. 

This  case  is  handled  in  step  4,  below. 

3  -  If  the  (originally)  calculated  value  of  F  is  negative, 

’  J  out , X 

=  0  and 

out,x 

the  resulting  system  is  solved  for  F 


out,x 

then  the  x  diamond  relation  is  replaced  by  the  fixup  F 


out,y 


„  Q  +  F.  a  +  F.  (b  -  1/2) 

F  =  m,x  m,y 

out,y  - * - 


(V-31) 


b  +  1/2 


If  F  is  not  negative,  then  compute  F  from  the  y  diamond  relation, 

equation  (V-27),  which  completes  the  calculations.  Otherwise,  both 

outgoing  fluxes  require  fixups.  This  case  is  handled  in  step  4,  next. 

4  -  Since  both  outgoing  fluxes  have  been  negative,  despite 

application  of  a  fixup  of  one  of  them,  the  final  alternative  is  to 

replace  both  diamond  relations  by  fixup  constraints:  F  =0  and 

out,x 

Fout  y  =  0  an^  compute  Fc  directly  from  the  balance  equation: 
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Q  + 


F. 
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F.  b 
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V-l 


6.  Step  Characteristic 

The  step  characteristic  method  for  xy-geometry  was  developed 
by  Lathrop  and  compared  to  several  other  positive  schemes  [Ref.  7]. 
Figure  V-l  illustrates  a  single  space  mesh  unit  cell,  of  dimensions  Ax 
by  Ay,  not  necessarily  square.  The  streaming  direction  is  assumed  to 
lie  in  the  principal  octant,  so  that  particles  enter  the  cell  from  the 
left  and  bottom  and  exit  on  the  right  and  top. 


Rather  than  assume  the  auxiliary  relations  directly,  as  in  the 


previous  methods,  the  step  characteristic  method  uses  analytic 

A 

integration,  in  the  fixed  single  direction  0  ,  to  find  the  flux 
tnroughout  the  cell.  The  assumptions  are: 

1  -  the  source  term,  Q,  is  constant  throughout  the  cell 

2  -  the  flux  entering  on  the  left,  F.  ,  is  a  constant  (with 

m,x 

respect  to  y)  along  the  left  face  of  the  cell 

3  -  the  flux  entering  on  the  bottom,  F.  ,  is  a  constant 

m,y 

(with  respect  to  x)  along  the  bottom  face  of  the  cell. 

The  solution  for  the  flux  is  analytic  throughout  the  cell,  except 

that  it  is  discontinuous  across  the  corner  characteristic  line,  which 

extends  from  the  lower  left  corner  in  direction  <(>  across  the  cell 

(unless  F.  =  F.  ).  The  required  outputs  of  the  method  are  the  cell 
m,x  m,y 

average  flux,  F  ,  and  the  exiting  fluxes,  F  ^  and  F  .  The  exiting 
c  3  out,x  out,y 

fluxes  are  obtained  by  analytically  averaging  the  solution  for  the  flux 
along  the  right  and  top  faces.  This  provides  a  formula  for  the  outbound 
fluxes  in  terms  of  the  inbound  fluxes,  the  source  (Q) ,  and  the  optical 
thicknesses  of  the  cell  («and  g).  The  cell  average  flux,  F^,  could  be 
found  by  analytically  averaging  the  solution  over  the  cell,  but  it  is 
computationally  more  efficient  to  find  it  directly  from  the  balance 
equation.  The  two  approaches  are  algebraically  equivalent,  but  not 
necessarily  numerically  equivalent,  due  to  rounding  error.  Using  the 
balance  equation  ensures  that  the  method  is  numerically,  as  well  as 
conceptually,  conservative.  The  step  characteristic  method  is  a  variable 
weight  method,  in  that  there  are  different  sets  of  formulas  used 
depending  upon  whether  the  corner  characteristic  line  intersects  the 
right  face,  the  top  face,  or  the  top  right  corner  of  the  cell. 
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If  the  corner  characteristic  line  intersects  the  right  face  of  the 

cell,  as  shown  in  figure  V-l,  then  <J>  <  <j>  . .  =  arctan  (Ay/A  x) .  An 

diag 

equivalent  condition  is  that  a  <8  .  In  this  case,  the  right  cell-edge 
flux  depends  on  both  inbound  fluxes,  but  the  top  flux  does  not  depend  on 
the  bottom  flux.  The  formulas  are 


F  =  Q  +  (1-  ab)  e"“(F.  -Q)  +  (l-e~a)  b  (F.  -Q) 

out,x  m,x  *  in,y 

F  =  Q  +  (l-e_0t)  a  (F,  -Q) 

out,y  m,x 


(V-33) 

(V-34) 


where  a  =  1/a  and  b  =  1/8  .  If  the  corner  characteristic  line 
intersects  the  top  right  corner,  i.e.  if  a  =8  ,  then  the  top  flux 
depends  only  on  the  left,  and  the  right  depends  only  on  the  bottom.  The 
formulas  are 


F 

out,x 

F 

out,y 


Q  +  (l-e_B)  b  (F  -Q) 
in  t  y 

Q  +  (l-e-a)  a  (F.n>x-Q) 


(V-35) 
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The  remaining  case  is  that  a  >6,  so  that  the  formulas  are 


F 

out,x 

F  . 
out,y 


-8 

Q  +  (1-e  )  b  (F,  -Q) 

m,y 

-8  »  -8 
Q  +  (1-e  )  a  (F,  -Q)  +  (l-8a)  e  (F.  -Q) 

in,x  m,y 


(V-37) 
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In  each  case,  the  cell  average  flux  is  given  by 


The  step  characteristic  spatial  quadrature  is  advantageous  for  use 
in  the  discrete  elements  method  for  both  the  auxiliary  and  main  fluxes. 
It  is  unconditionally  stable  and  convergent,  as  well  as  being  a  positive 
method.  These  features  provide  the  smooth  results  needed  for  the 
auxiliary  fluxes,  but  with  much  better  accuracy  than  the  step 
quadrature.  The  step  characteristic  quadrature  also  has  the  accuracy 
needed  for  the  main  fluxes.  However,  the  accuracy  of  the  method  depends 
upon  the  relation  between  the  streaming  azimuth  and  the  cell  shape.  If 
the  streaming  direction  is  at  or  near  the  cell  diagonal,  the  method  has 
excellent  accuracy,  but  it  provides  comparatively  poor  accuracy  for 
directions  which  are  far  from  the  diagonal.  The  nature  of  these 
inaccuracies  is  considered  in  detail  in  the  next  chapter. 

C.  Summary 

This  chapter  has  reviewed  the  spatial  quadrature  methods  used  in 
one-  and  two-dimensional  discrete  ordinates  and  for  evaluation  of  the 
discrete  elements  method.  The  step  method  is  computationally 
inexpensive,  and  is  inaccurate  but  positive  and  smooth.  The  diamond 
difference  method  is  computationally  moderately  expensive  and  is 
accurate  in  computing  global  quantities,  such  as  total  absorption  rate, 
but  is  less  accurate  for  pointwise  quantities,  such  as  the  scalar  flux 
as  a  function  of  position.  This  inaccuracy  is  because  diamond  difference 
is  not  a  positive  method,  and  hence  predicts  negative  fluxes  (if  the 
space  mesh  is  too  coarse)  or  spatially  oscillating  fluxes,  in  any  case. 
The  diamond  difference  method  may  be  used  with  negative  flux  fixups, 
which  introduce  nonlinearity  into  the  method,  and  reduce  accuracy  of 


global  quantities,  but  which  prevent  negative  fluxes  and  so  improve  the 


calculation  of  pointwise  quantities.  The  step  characteristic  method  is 


computationally  relatively  expensive  because  of  the  exponential  function 
it  uses,  but  provides  both  smoothness  and  accuracy.  In  one  dimension, 
the  linear  characteristic  method  is  also  selected  for  testing.  It  is 
higher  order  than  the  step  characteristic  and  has  improved  accuracy, 
especially  with  relatively  coarse  spatial  meshes. 

Later  chapters  report  the  results  of  testing  the  discrete  element 
method  with  various  combinations  of  methods  of  spatial  quadrature  used 
for  the  auxiliary  and  main  fluxes.  In  order  to  interpret  these  results, 
it  is  necessary  to  consider  the  characteristics  of  the  spatial 
quadratures  as  well  as  those  of  the  angular  quadratures.  This  is  the 
reason  for  the  detailed  review  presented  here.  The  interaction  of  the 
spatial  and  angular  quadratures  is  similarly  important,  and  is 
considered  in  the  next  chapter. 
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VI .  Interaction  of  Angular  and  Spatial  Quadrature: 


Ray  Effects  and  Quasi-Ray  Effects 

The  objective  of  this  research  is  to  develop  and  evaluate  a  method 
of  numerical  neutron  transport  which  handles  difficult  problems  better 
than  discrete  ordinates  and  which  reduces  ray  effect.  The  discrete 
elements  method  attempts  this  by  coupling  the  angular  and  spatial 
quadrature  and  by  using  a  composite  angular  quadrature.  In  order  to 
accurately  interpret  the  results  of  numerical  testing,  it  is  first 
necessary  to  consider  what  sorts  of  errors  are  attributable  to  the 
angular  quadrature  in  SN  and  what  sorts  are  due  to  the  spatial 
quadrature,  and  how  the  two  interact.  The  purpose  of  this  chapter  is  to 
investigate  these  deficiencies  and  interactions.  The  approach  is  to 
solve  a  simple  problem,  a  square  source  in  a  square  non-scattering 
absorber,  using  combinations  of  analytic  and  discrete  angular  and 
spatial  quadrature. 

A.  T^e  Square -in-a-Square  Problem:  Without  Scatter 

Figure  VI-1  shows  the  test  problem  and  a  benchmark  solution.  The 

problem  is  a  4  cm  by  4  cm  square  absorber  with  an  absorption  cross- 

section  of  0.75  cm  *  and  scatter  cross-section  of  zero,  located  in  a 

2 

vacuum.  A  source  of  strength  1  n/cm  /sec  is  distributed  throughout  the 
central  1  cm  by  1  cm  subregion.  Only  the  upper  right  quadrant  of  the 
problem  is  actually  solved.  The  rest  of  the  problem  is  represented  by 
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reflecting  boundaries.  Also  shown  is  a  benchmark  solution  for  the  y 


component  of  scalar  current,  J^,  along  the  top  edge  of  the  square.  This 
solution  used  Monte  Carlo  spatial  and  azimuthal  integration  and  2-point 
Gauss-Christof fel  polar  integration  of  a  Green's  function  kernel.  The 
curve  is  obtained  from  16  data  points,  each  representing  the  average 
current  in  an  interval  1/8  cm  in  width. 

This  benchmark  solution  shows  several  qualitatively  expected 
features.  The  current  is  nowhere  zero.  The  current  is  monotone 
decreasing,  as  a  function  of  x  (along  y  =  2  cm).  The  current  should  have 
zero  slope  with  respect  to  x  at  x=0  (by  symmetry),  and  the  solution 
approximates  this  feature.  The  slight  lack  of  smoothness  is  due  to  the 
variance  of  this  Monte  Carlo  solution,  which  used  200,000  particles. 


Jy(x,y«2) 
(Graphed 
at  right) 


Fig.  VI-1:  Square-in-a-Square ,  Non-scattering  Absorber  Problem 


B.  Angular  Truncation  Error  and  Ray  Effect 

Any  finite  numerical  representation  of  the  angular  dependence  of 
the  flux  will  suffer  truncation  error.  This  results  in  inaccuracy  of  the 
computed  scalar  fluxes,  which  feeds  back  through  the  scattering  source 
in  the  transport  equation  to  cause  inaccuracy  in  the  spatial  variation 
as  well.  One  of  the  objectives  of  the  discrete  element  representation  is 
to  reduce  this  error  and  thus  improve  the  accuracy  of  computed 
quantities . 

Ray  effect,  however,  is  a  particular  sort  of  systematic  error  which 
causes  computed  results  to  be  qualitatively  wrong,  as  well  as 
quantitatively  inaccurate.  The  term  "ray  effect"  is  often  used  loosely 
to  describe  any  such  qualitative  deficiencies.  For  the  purposes  of  this 
discussion,  a  more  precise  definition  is  needed. 

Ray  Effect:  The  presence  of  qualitatively  unreasonable 
numerical  results  caused  by  the  use  of  a  discrete 
angular  representation  of  the  directional  flux,  which 
would  be  present  even  if  the  spatial  quadrature  were 
performed  analytically,  or  with  a  vanishingly  fine  mesh. 

As  an  example  of  ray  effect,  figure  VI-2  shows  results  for  the  test 
case  of  figure  VI-1.  These  are  Green's  function  solutions  with  discrete 
ordinates  angular  quadrature,  evaluated  at  the  midpoints  of  16  intervals 
of  width  1/8  cm  along  the  top  edge  of  the  problem.  The  S  and  S 
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solutions  show  the  results  of  projecting  a  single  ray  from  the  source 

region  to  the  edge  of  the  absorber.  These  solutions  are  neither  non-zero 

everywhere  nor  monotone  decreasing.  However,  comparison  with  the 

benchmark  solution  shows  to  be  essentially  converged  in 
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azimuthal  quadrature. 
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C.  Quasi-Ray  Effect  and  Spatial  Truncation  Error 

As  with  the  angular  representation,  the  finite  numerical 
representation  of  the  spatial  variation  of  the  fluxes  results  in 
truncation  error  causing  quantitative  inaccuracy  of  computed  results. 
Although  less  generally  recognized,  systematic  errors  can  also  occur, 
resulting  in  qualitative  inaccuracy.  These  inaccuracies  are,  in 
practice,  quite  difficult  to  distinguish  from  those  caused  by  the 
angular  discretization,  but  the  conceptual  distinction  is  important.  In 
view  of  this,  the  term  quasi-ray  effect  is  defined  for  use  here. 

Quasi-Ray  Effect:  The  presence  of  qualitatively  unreasonable 
numerical  results  caused  by  the  use  of  a  discrete 
spatial  representation  of  the  directional  flux,  which 
would  be  present  even  if  the  angular  quadrature  were 
performed  analytically,  or  with  a  vanishingly  fine  mesh. 

As  an  example  of  quasi-ray  effect,  the  test  case  of  figure  VI-1  is 
solved  by  S__  using  discrete  spatial  quadratures.  Since  the  angular 
quadrature  is  essentially  converged,  any  appearance  of  ray  effect  should 
be  attributable  to  the  spatial  quadrature.  The  spatial  mesh  used  is  16 
by  16  cells  of  size  1/8  cm  by  1/8  cm.  The  resulting  curves  of  are 
thus  comparable  to  those  presented  above.  Figure  VI-3A  shows  the  results 
for  two  spatial  quadrature  methods,  diamond  difference  with  negative 
flux  fixup  (DDF)  and  step  characteristic  (SC).  With  a  total  cross- 
section  of  0.75  cm  \  the  cells  are  slightly  less  than  1/10  mean  free 
path  in  height  and  width,  and  this  is  a  fine  enough  mesh  that  the  DDF 
curve  nearly  matches  the  cori esponding  analytic  curve  in  figure  VI-2. 
The  spatial  oscillations  of  DDF,  like  noise,  have  averaged  out. 
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The  SC  solution  in  figure  V1-3A  shows  a  distinct  quasi-ray  effect. 

The  flux  directly  above  the  source  region  is  too  high  and  that  just 

adjacent  is  too  low.  This  comes  about  in  the  following  way.  Consider  one 

of  the  azimuthal  quadrature  directions  in  the  ,  „„  set  which 

represents  streaming  of  particles  upward  and  to  the  left,  as  shown  in 

figure  VI-3B.  If  the  SC  method  were  to  accurately  represent  a  collimated 

beam  of  particles  streaming  in  this  direction,  starting  with  a  source 

only  in  cell  7,  then  the  flux  would  step  through  the  space  mesh 

obliquely,  as  does  a  knight  on  a  chess  board,  moving  from  cell  7  through 

cell  6  of  that  figure,  so  that  none  of  the  flux  would  reach  cells  3  and 

9.  However,  because  the  SC  method  assumes  a  constant  flux  distribution 

along  each  cell  interface,  part  of  the  flux  is  "averaged"  to  the  right 

along  the  7-8  interface,  and  again  along  the  8-9  interface,  arriving, 

incorrectly,  in  cell  9.  Similarly,  some  of  the  flux  crosses  through  cell 

4  into  cell  5,  and  is  "averaged"  to  the  left  along  the  4-5  interface, 

thence  crossing  through  cell  2  into  cell  3.  Thus,  some  of  the  flux  goes 

off  to  the  left  like  a  bishop,  while  most  of  the  flux  goes  up  the  grid 

like  a  rook.  Rather  than  a  collimated  beam,  the  streaming  is  modelled  as 

a  broad  swath  which  more-or-less  averages  out  to  motion  in  the  intended 

direction.  The  closer  that  direction  is  to  diagonal,  the  better 

collimated  the  beam;  but  for  directions  far  from  the  diagonal,  the  flux 

moves  straight  up  the  grid  with  a  small  smear  out  to  one  side.  Since  the 

S.  quadrature  includes  many  such  oblique  angles,  this  error 
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accumulates  systematically,  overestimating  the  flux  directly  above  the 
source.  The  flux  is  decreased  just  outside  this  band,  i.e.  at  x  just 
greater  than  1/2  cm,  since  the  particles  that  should  have  arrived  there 
have  increased  the  flux  above  the  source,  instead.  Similarly,  transport 
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along  the  cell  diagonals  (like  a  bishop)  accumulates  so  that  the  flux 
diagonally  opposite  the  source  (near  x  =  2  cm)  is  also  overestimated. 

The  overall  result  is  qualitative  error  which  is  similar  to  what 
one  would  expect  for  ray  effect  if  a  quadrature  set  which  included  only 
the  verticle  and  diagonal  directions  were  used.  This  is  the  quasi-ray 
effect  defined  above. 

D.  Combined  Angular  and  Spatial  Truncation  Errors 

In  practice,  both  the  spatial  and  angular  quadratures  are  discrete 
so  that  it  is  difficult  to  distinguish  ray  effect  and  quasi-ray  effect. 
The  use  of  analytic  spatial  or  angular  quadrature  is  normally 
infeasible.  Only  the  simple  example  chosen  here,  without  scatter,  has 
made  it  possible  to  demonstrate  the  independent  nature  of  each. 

Figure  VI-4A  shows  the  S  solutions  for  step,  DDF,  and  SC 
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spatial  quadratures  for  comparison  with  the  analytic  spatial  quadrature 
of  figure  VI-2.  The  step  method  results  are  very  smooth  and  show  the 
presence  of  flux  throughout  the  absorber  region,  for  reasons  similar  to 
those  discussed  above  with  reference  to  the  SC  method.  The  DDF  results 
resemble  the  analytic  solution  well  in  that  regard  (no  current  for  x  <  1 
cm),  but  show  the  lack  of  smoothness  typical  of  the  method.  The  SC 
results  show  excellent  accuracy  (of  spatial  quadrature)  for  this  case  of 
diagonal  streaming. 

Figure  VI-4B  shows  the  S  solutions  for  DDF  and  SC  spatial 
quadratures.  Compared  with  figure  VI-2,  the  DDF  results  are  generally 
correct,  but  lack  smoothness.  The  SC  results  are  smooth,  but  show  the 
beam  spreading  described  above  for  this  non-diagonal  streaming 


direction. 


Carlson  and  Lathrop  have  remarked  upon  the  interaction  of  spatial 


and  angular  quadrature. 

"Experience  has  shown  that  errors  involved  in  spatial 
and  angular  quadrature  are  interdependent.  Qualitatively,  the 
error  surface  is  like  a  valley  between  two  ridges.  If  error 
is  plotted  against  order  of  angular  quadrature  along  one  axis 
and  order  of  spatial  quadrature  along  an  orthogonal  axis,  the 
ridges  of  the  surface  lie  above  these  axes.  Hence,  if  a 
calculation  gives  a  result  in  the  error  valley,  both 
quadratures  must  be  refined  to  remain  in  this  error  valley. 

On  the  other  hand,  for  a  given  spatial  mesh,  refining  the 
angular  quadrature  may  actually  increase  the  error,  and 
conversely."  [Ref.  4:35] 

The  examples  presented  in  this  chapter  tend  to  explain  these 
observations.  A  low  order  angular  quadrature  should  suffer  from  strong 
ray  effect,  due  to  the  use  of  a  few  narrow  beams  of  streaming  particles; 
but  the  spatial  quadrature,  for  a  similarly  coarse  mesh,  broadens  out 
these  narrow  beams,  reducing  the  apparent  ray  effect.  Arbitrarily 
refining  only  the  spatial  mesh  would  narrow  down  the  beams,  letting  the 
ray  effect  be  seen.  Conversely,  a  coarse  spatial  mesh  should  suffer  from 
strong  quasi-ray  effect,  but  a  correspondingly  coarse  angular  mesh 
prevents  this.  With  only  a  few  angles,  the  spatial  errors  seem  random, 
and  are  considered  as  mere  truncation  error.  But  arbitrarily  refining 
only  the  angular  mesh  causes  accumulation  of  this  systematic  error, 
letting  the  quasi-ray  effect  be  seen. 

E.  Summary 

This  chapter  has  demonstrated  some  characteristics  of  the  discrete 
ordinates  method. 

1  -  Spatial  quadrature  schemes  smear  out  the  angular  rays 

2  -  Use  of  a  limited  number  of  angular  directions  often  ameliorates 


the  quasi-rays  of  the  spatial  scheme 


Success  of  S„  may  be  due,  in  part,  to  the  combination  of  these  two 
N 

effects.  On  the  other  hand,  the  general  truncation  error  implicit  in  the 
method  has  precluded  application  to  some  types  of  problem,  such  as 
vacuum  ducts. 

The  discrete  elements  method  avoids  ray  effects  not  by  smearing  out 

fixed  rays,  but  rather  by  steering  the  rays.  This  means  of  coupling  the 

spatial  and  angular  representations  then  requires  a  spatial  quadrature 

scheme  which  will  accurately  propagate  the  element  flux  as  a  collimated 

beam  in  the  steered  direction.  Recent  progress  has  been  made  in 

developing  higher  order  spatial  quadrature  schemes,  such  as  the  linear 

characteristic  and  linear  nodal  methods,  which  should  be  of  value  in 

this  regard.  The  discrete  elements  method  may  also  be  viewed  as  a  higher 

order  angular  quadrature  scheme,  which  should  complement  these  new 

spatial  quadratures.  Success  of  the  discrete  elements  method  should  be 

judged  both  on  its  ability  to  produce  the  results  of  a  high  order 

discrete  ordinates  calculation  with  the  use  of  fewer  rays,  and  on  its 

ability  to  handle  difficult  problems  more  accurately  than  S  .  The 

N 

following  chapters  present  the  results  of  numerical  testing  of  this 
ability. 
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VII.  Versions  of  the  Discrete  Elements  Method 


The  previous  chapters  have  shown  that  there  are  several  degrees  of 
freedom  in  designing  a  discrete  elements  scheme:  quadrature  rule  for  the 
individual  elements,  angular  mesh  type  and  size,  spatial  quadrature  for 
the  auxiliary  and  main  flux  calculations  (not  necessarily  the  same),  and 
the  coupling  between  the  auxiliary  and  main  fluxes. 

A.  Specification  of  a  Discrete  Elements  Numerical  Scheme 

Notation  to  identify  the  particular  combination  of  choices 

constituting  a  specific  discrete  element  method  is  defined  as  follows: 

Prefix  -  the  general  notation,  L„,  is  preceded  by  an  abbreviation 
-  N 

indicating  the  quadrature  rule  used  within  each  discrete  element: 

SR  -  Simpson's  rule 

G3  -  Gauss-Legendre  three-point  quadrature 
NC  -  Newton-Cotes  (five-point)  rule 
Subscripts  -  the  general  subscript,  N,  is  replaced  by  a  subscript 
indicating  the  type  of  angular  mesh  and  number  of  mesh  intervals.  These 
subscript  conventions  were  introduced  in  chapter  III. 

Superscripts  -  the  superscripts  indicate  the  spatial  quadrature 
method(s)  used  for  the  auxiliary  and  main  fluxes.  The  abbreviations  for 
the  superscripts  are: 

Step  -  step  method 
DD  -  diamond  difference 

DDF  -  diamond  difference  with  negative  flux  fixup 
SC  -  step  characteristic 
LC  -  linear  characteristic 


SC /DDF 

Thus,  a  notation  of  G3-L_„  would  indicate  a  discrete 
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elements  scheme  with  three-point  Gauss-Legendre  quadrature  in  azimuth 
within  each  element  of  a  three  cell  equal  weight  azimuthal  mesh  and  with 
two-point  Gauss-Christof fel  polar  quadrature,  where  the  auxiliary  fluxes 
are  calculated  by  the  step  characteristic  method  and  the  main  fluxes  by 
diamond  difference  with  negative  flux  fixup.  Portions  of  the  notation 
which  are  determined  by  the  context  of  a  discussion  are  usually  omitted, 
for  simplicity. 

B.  Coupling  of  Auxiliary  and  Main  Fluxes 

In  computing  the  auxiliary  fluxes  by  the  discrete  ordinates  method, 
it  is  possible  to  either  use  the  scalar  flux  from  the  main  calculation 
or  to  compute  a  separate  scalar  flux  directly  from  the  auxiliary  fluxes. 
In  the  former  case,  the  auxiliary  fluxes  are  coupled  to  the  main  fluxes 
by  feedback  through  the  source  term.  In  the  latter  case,  the  auxiliary 
fluxes  are  uncoupled  from  the  main  fluxes.  (Note  that  the  main  fluxes 
are  always  coupled  to  the  auxiliary  fluxes  through  the  estimation  of  the 
streaming  directions  of  the  main  fluxes,  which  is  the  object  of  the 
exercise.)  Coupling  would  seem  to  offer  more  accurate  auxiliary  fluxes, 
and  hence  more  accurate  steering  of  the  main  fluxes,  assuming  stability 
is  satisfactory.  On  the  other  hand,  uncoupling  preserves  the  linearity 
of  the  auxiliary  flux  calculation,  ensuring  convergence.  Uncoupling  also 
means  that  the  streaming  directions  of  the  main  fluxes  are  fixed,  in  the 
sense  of  being  determined  independently,  so  that  the  convergence  of  the 
main  calculations  is  also  assured. 
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The  relative  merits  of  coupling  vs.  uncoupling  were  evaluated  in 
one-dimensional  L„  and  coupling  was  consistently  the  more  accurate 
technique.  Uncoupling  speeded  convergence  only  slightly,  for  the 
problems  considered.  Uncoupling  could  allow  the  streaming  directions  to 
be  computed  in  advance  and  stored  for  use  in  the  main  calculations 
(although  this  would  require  rather  vast  storage  resources  for  practical 
problems),  but  the  decrease  in  computational  effort  of  this  scheme 
proved  negligible.  It  was  concluded  that  coupling  is  the  scheme  to  use. 
If  poor  convergence  is  encountered,  convergence  can  be  improved  by 
refining  the  spatial  mesh  or  using  a  smoother  spatial  method,  rather 
than  by  uncoupling.  The  results  and  conclusions  presented  in  this 
dissertation  refer  implicitly  to  the  coupled  form  of  discrete  elements. 

C.  Optimization 

For  computer  implementation,  a  numerical  method  may  be  optimized  in 
two  important  ways:  with  respect  to  execution  time  and/or  with  respect 
to  storage.  The  dollar  cost  of  a  calculation  may  be  dominated  by  either 
of  these  factors,  depending  on  the  computer  facility.  One  of  the 
objectives  of  this  research  is  to  identify  the  "best"  discrete  element 
method,  based  on  accuracy  vs.  execution  time  with  near  minimum  storage. 
The  next  chapter  presents  execution  time  and  storage  requirements  for 
two-dimensional  versions  of  the  discrete  elements  method. 
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This  chapter  compares  the  computer  execution  time  and  memory 


requirements  of  the  discrete  elements  and  discrete  ordinates  methods. 
This  information  facilitates  cost-effectiveness  comparisons  of  the  two 
methods  and  establishes  the  feasibility  of  the  discrete  elements  method. 

A.  Problem  and  Quadrature  Dimensions 

The  following  memory  requirements  apply  to  a  problem  mesh 
consisting  of  I  cells  in  the  x  direction  by  J  cells  in  the  y  direction, 
and  K  polar  directions  or  elements  by  L  azimuthal  directions  or 
elements.  The  IJ  space  cells  are  partitioned  into  R  regions  of  constant 
material  and  source  parameters,  not  necessarily  contiguously  located. 

B.  Discrete  Ordinates  Memory  Requirements 

The  discrete  ordinates  method  requires  the  following  arrays: 

Spatial  Arrays  -  Three  arrays  of  size  IJ  are  required  for  the  all¬ 
angle  flux  for  the  previous  iteration,  F  ,  ,,  the  scalar  flux  for  the 

old 

iteration  in  progress,  F  w,  and  the  indexing  array  which  records  to 
which  region  each  cell  belongs.  Three  arrays  of  size  R  are  required  for 
the  material  constants,  c,  o,  and  Qq  (=  S/o).  One  cell  spacing  array  of 
size  I  is  required  for  Ax,  and  one  of  size  J  for  Ay,  assuming  non- 
uniform  cell  sizes. 

Angular  Arrays  -  For  minimum  storage,  two  arrays  of  size  K  are 
required:  one  for  sin(0)  and  one  for  the  polar  weights,  W^.  Two  arrays 
of  size  L  are  required  for  the  azimuthal  angles  and  weights.  Execution 
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speed  is  slightly  enhanced  if  arrays  of  size  KL  are  used  for  y,  n,  and 


Cell  Boundary  Flux  Arrays  -  If  two  vacuum  boundaries  are  available, 
then  the  boundary  values  of  the  directional  fluxes  need  not  be  stored 
between  iterations.  Within  an  iteration,  two  arrays  of  size  I  and  two  of 
size  J  are  required  for  the  intercell  boundary  fluxes.  If  there  are  no 
vacuum  boundaries,  then  KL  sets  of  these  arrays  are  required  to  save  the 
needed  boundary  values  of  the  directional  fluxes  between  iterations. 

Optical  Thickness  Arrays  -  Although  not  required,  execution  speed 
is  improved  if  the  optical  thicknesses  of  cells  in  each  region  and  in 
each  quadrature  direction  are  stored.  Assuming  a  uniform  space  mesh 
within  each  region,  this  requires  two  arrays  of  size  RKL,  otherwise  it 
would  require  two  arrays  of  size  IJKL.  This  is  probably  impractical 
except  in  the  former  case  and  with  relatively  small  R,  and  it  is  assumed 
that  these  arrays  are  not  used. 

Minimum  Memory  Requirement: 

With  Vacuum  Boundaries:  31 J  +  3R  +  3(I+J)  +  2(K+L) 

Without  Vacuum  Boundaries:  31 J  +  3R  +  (2KL+1) (I+J)  +  2(K+L) 

C.  Discrete  Elements  Memory  Requirements 

The  discrete  elements  method  requires  the  storage  described  above 
for  the  discrete  ordinates  method,  with  the  element  main  fluxes  taking 
the  place  of  the  ordinate  fluxes,  and  has  additional  storage 
requirements  for  the  auxiliary  fluxes.  Assuming  the  element  quadrature 
rule  has  A  auxiliary  fluxes  per  element,  the  following  minimum  of  arrays 
are  required: 

With  two  adjacent  vacuum  boundaries,  2A  arrays  of  size  I  and  of 
size  J  are  needed  for  the  intercell  auxiliary  fluxes  and  2A  arrays  of 
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size  K  and  of  size  L  are  needed  for  the  auxiliary  direction  quadrature 
parameters.  Without  vacuum  boundaries,  KL  sets  of  2A  arrays  of  size  I 


and  of  size  J  are  needed  for  the  auxiliary  flux  boundary  values,  in 
addition  to  the  auxiliary  quadrature  parameter  arrays. 

Minimum  Memory  Requirement: 

With  Vacuum  Boundaries: 

3IJ  +  3R  +  (2A+3) (I+J)  +  (2A+2) (K+L) 

Without  Vacuum  Boundaries: 

31 J  +  3R  +  [2KL(A+1)+1]  (I+J)  +  2(A+1)(K+L) 


Fig.  VIII-1:  Quadrature  Set  for  SR-L„_ 

-  2G ,  3 

Additional  memory  may  be  used  to  speed  execution  with  Simpson's 

rule  or  Newton-Cotes  rule,  as  mentioned  in  chapter  III.  As  an  example, 

figure  VIII-1  shows  the  quadrature  set  for  the  SR-L  scheme  for  the 

2G,  3 

principal  octant.  The  circles  mark  the  element  center  auxiliary 
directions;  the  squares  mark  the  element  edge  auxiliary  directions  at 
the  edges  of  the  octant;  the  diamonds  mark  the  element  edge  auxiliary 
directions  internal  to  the  octant.  Each  of  these  last  four  directions  is 
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common  to  two  elements,  and  could  be  computed  (throughout  the  space 


mesh)  only  once,  with  the  first  element  walk,  and  reused  for  the  second 

element  walk.  This  results  in  a  10%-15%  savings  in  execution  time  for 

this  quadrature.  The  cost  in  storage  is  an  additional  array  of  size-  IJ, 

if  no  vacuum  boundaries  are  used,  and  two  such  arrays  if  vacuum 

boundaries  are  used  (since  octants  are  done  in  pairs:  a  single  elements 

is  walked  in  from  vacuum  boundary  to  non-vacuum  boundary,  reflected,  and 

walked  back  out  to  the  vacuum  boundary).  To  take  full  advantage  of  this 

quadrature  set  overlap  in  SR-L  requires  (K+l)  sets  of  these  arrays.  K 

K,  L 

sets  are  required  for  the  directions  marked  with  triangles  in  figure 

VIII-2,  an  example  showing  SR-L  ,  and  one  set  is  needed  for  the 

2/3 

diamonds.  This  provides  about  a  15%  savings  for  this  example,  and  up  to 
about  25%  for  a  fine  angular  mesh. 
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Fig.  VIII-2:  Quadrature  Set  for  SR-L„  „ 

2  /  J 

This  memory  requirement,  (K+l)  IJ,  may  be  prohibitive,  but  it  can  be 

quite  reasonable  when  compared  to  the  GIJ  storage  (for  F  ,  )  required 

old 

for  G  energy  groups,  since  G  is  often  large  and  K=4  is  probably  maximum. 
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D.  Example  Comparison  of  Storage  Requirements 

Table  VIII-1  shows  the  memory  requirements  for  an  example  case  of  a 

100  cell  by  100  cell  space  mesh  with  25  regions.  The  minimum  memory 

requirements,  measured  in  array  elements,  are  given  for  two  situations, 

with  and  without  vacuum  boundaries.  The  formulas  for  conventional  are 

N 

slightly  different  than  those  above,  and  the  result  is  shown  for  S^. 


Method 

With  2  Vacuum  Boundaries 

No  Vacuum  Boundaries 

S16 

30,783 

44,783 

S2G,3R  30,685  32,685 


L2G,3 

31,915 

39,915 

L2,3 

32,735 

44,735 

Table  VIII-1:  Example  Case  Memory  Requirements 

With  two  adjacent  vacuum  boundaries,  the  extra  memory  required  for 
the  discrete  elements  method  is  negligible.  Without  vacuum  boundaries, 
it  is  larger,  but  not  problematical.  For  very  large  spatial  mesh 
problems,  where  memory  may  be  a  limiting  factor,  the  three  arrays  of 
size  IJ  dominate  all  the  schemes,  so  that  the  discrete  element  method  is 
essentially  equivalent  to  discrete  ordinates  in  storage  requirements. 

E.  Execution  Times  for  SN  Spatial  Quadratures 

The  execution  speed  of  the  discrete  ordinates  method  is  determined 
almost  entirely,  by  the  spatial  quadrature  subroutine.  Table  VIII-2  gives 
the  sopri  x  '  ..ate  execution  costs  of  elementary  mathematical  floating 
point  operations  and  function  evaluations  in  units  of  Cray-1  clock 
cycles.  A  short  accuracy  square  root  routine  of  approximately  six 
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decimal  place  accuracy  is  estimated  to  be  achievable  using  65  clock 
cycles.  A  short  accuracy  (four  to  six  decimal  places)  routine  for 
evaluating  both  the  sine  and  cosine  costs  46  cycles  using  the  following 
formulas  (Ref.  1:76): 

v  =  u  *  u  (VIII-1) 
cos(u)  =  (.03705  *  v  -  .49670)  *  v  +  1.0  (VIII-2) 
sin (u)  =  ((.00761  *  v  -  .16605)  *  v  +  1.0)  *  u  (VIII-3) 


Operation:  +  -  *  /if  exp  sqrt  sin  &  cos  (both) 

Clock  Cycles:  4  4  5  20  8  128  65  46 

Table  VIII-2:  Cray-1  Floating  Point  Execution  Times 

Estimates  of  the  execution  costs  of  the  two  dimensional  spatial 
quadrature  schemes  of  chapter  V,  based  on  the  above  clock  times,  are 
presented  in  table  VIII-3.  Larsen  and  Alcouffe  [Ref.  5]  have  recently 
developed  a  linear  characteristic  method  for  xy  geometry  with  a  cost  of 
about  600  clock  cycles.  The  diamond  difference  method  requires  more  than 
the  indicated  102  clock  cycles  for  cells  where  fixups  are  required;  the 
step  characteristic  method  requires  slightly  less  than  267  cycles  if  the 
quadrature  direction  lies  exactly  on  the  space  cell  diagonal.  The 
numbers  given  are  generally  representative,  however. 


XY  Spatial  Quadrature: 

Step 

DDF 

SC 

LC 

Clock  Cycles  per  Call: 

56 

102 

267 

600 

Table  VIII-3:  Discrete  Ordinates  Spatial  Quadrature  Subroutine  Costs 


F.  Execution  Times  for  Streaming  Direction  Averaging 

The  cost  of  execution  of  a  discrete  elements  scheme  is  determined 
by  the  sum  of  three  sets  of  subroutine  calls: 

1  -  Auxiliary  flux  spatial  quadrature  (per  table  V1II-3) 

2  -  Streaming  direction  flux-weighted  averaging  (per  table  VIII-4) 

3  -  Main  flux  spatial  quadrature  (per  table  VIII-3) 

The  spatial  quadratures  are  computed  by  the  same  subroutine (s)  as  in 
discrete  ordinates,  although  different  schemes  can  be  used  for  the  main 
vs.  auxiliary  fluxes.  The  extra  cost  of  the  discrete  elements  method 
(besides  the  auxiliary  flux  calculations)  is  for  the  computation  of  the 
flux-weighted  streaming  directions.  These  costs  are  presented  in 
table  VIII-4,  below,  and  assume  the  use  of  short  accuracy  square  root, 
sine  and  cosine  algorithms.  Values  are  given  for  those  xy  geometry 
methods  which  were  tested.  An  alternative  scheme  of  directly  averaging 
the  x-  and  y-direction  cosines  was  tested  and  abandoned.  It  avoided  the 
sine  and  cosine  needed  to  convert  the  averaged  (x,(j>)  to  (y,n),  but  cost 
167  clock  cycles  and  gave  unsatisfactory  numerical  results.  The  short 
accuracy  sine  and  cosine  fits  proved  to  have  no  degrading  effect  on  the 
overall  accuracy  of  the  discrete  elements  method. 


1 


Method: 

Simpson's  Rule  or  Gauss-Legendre  (3-point) : 

Azimuthal  averaging  only  (L  ) 

KG  f  L 

Azimuthal  and  polar  averaging  (L 
Newton-Cotes  (5-point)  Rule: 


K,L 


) 


Clock  Cycles: 


160 

288 


Azimuthal  averaging  only  (L 


186 


As  an  example  of  the  use  of  the  above  information,  consider  a 
SC  sc  sc 

comparison  of  S,^  and  G3-L__  ,  '  .  The  discrete  ordinates  scheme, 

lb  ZG t  6 

S  ,  uses  144  total  quadrature  directions,  for  a  cost  of  38,448  cycles 

per  space  cell.  The  G3-L  discrete  elements  method  uses  three  auxiliary 

fluxes,  one  direction  averaging,  and  one  main  flux  for  a  cost  of  1228 

cycles  per  element  per  space  cell.  With  24  elements  total,  the  cost  is 

29,472  cycles  per  space  cell.  Assuming  both  converge  equally  rapidly, 

SC  sc 

these  figures  give  the  relative  cost  of  the  two  methods:  G3-L  ' 

2G  r  3 

SC 

costs  about  77%  as  much  as  S,, 

16 

H.  Computer  Programs  for  Testing  the  Discrete  Elements  Method 

The  computer  programs  used  to  test  the  discrete  elements  method 
were  written  in  Fortran-IV  (with  extensions  similar  to  WATFIV)  and  run 
on  an  Intel-8088  based  microcomputer.  The  programs  were  specialized  to 
handle  the  test  cases  efficiently  and  so  are  not  of  general 
applicability.  For  this  reason,  and  because  the  programs  were 
straightforward  implementations  of  the  equations  and  algorithms 
described  in  this  report,  listings  of  the  program  code  are  not  included 
here,  but  are  on  file  at  the  Physics  Department  of  the  Air  Force 
Institute  of  Technology,  Wright-Patterson  AFB,  Ohio,  45433. 

The  discrete  ordinates  codes  were  validated  by  comparison  with 
sample  cases  run  for  that  purpose  by  E.  W.  Larsen  of  the  Los  Alamos 
National  Laboratory.  The  discrete  elements  codes  were  validated  by  use 
of  the  midpoint  rule  for  the  element  quadrature  and  comparison  with  the 
equivalent  product  quadrature  discrete  ordinates  calculations. 
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I.  Conclusions 


As  with  most  higher  order  methods  in  numerical  analysis, 
computer  implementation  of  the  discrete  elements  method  entails  three 
costs:  increased  program  complexity,  increased  storage  requirements,  and 
increased  execution  times.  The  remainder  of  this  report  will  deal  with 
the  question  of  whether  the  method  also  provides  increased  accuracy, 
and,  if  so,  increased  efficiency  on  the  basis  of  cost  vs.  accuracy. 

Chapter  II  presented  the  discrete  elements  algorithm,  and  showed 
that  the  essential  simplicity  of  structure  of  the  discrete  ordinates 
method  is  retained. 

This  chapter  has  shown  that  the  L,  method  has  only  modestly 
increased  storage  requirements,  if  a  vacuum  boundary  is  present  in  both 
the  x-  and  y-directions.  In  the  absence  of  vacuum  boundaries,  the 
storage  requirements  are  more  substantial,  but  still  about  equivalent  to 
SN  requirements.  The  storage  penalty  is  of  even  less  concern  in  multi¬ 
group  problems,  since  the  increase  only  applies  to  the  one  energy  group 
being  iterated  at  a  time. 

Execution  times  are  similarly  increased,  but  are  comparable  to 

discrete  ordinates  quadratures  currently  in  use,  such  as  S  .  Two 

lb 

approaches  to  minimizing  execution  costs  are  possible.  One  is  to  use  the 
least  expensive  of  spatial  quadrature  schemes  for  the  auxiliary  fluxes, 
since  most  of  the  cost  of  the  method  is  for  these  calculations.  The 
other  approach  is  to  use  the  highest  order  of  spatial  quadratures  for 
both  the  main  and  auxiliary  fluxes  so  that  the  spatial  mesh  may  be 
coarse.  The  next  chapter  evaluates  these  approaches  for  the  case  of  one¬ 
dimensional  (slab)  geometry,  since  high-order  spatial  quadratures  are 
readily  employed  in  this  geometry. 
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IX.  Test  Cases  in  One-Dimensional  Geometry  and  Results 


This  chapter  compares  the  performance  of  Lvl  and  S._  schemes  for  two 

N  N 

problems  in  one-dimensional  Cartesian  geometry.  Evaluations  in  xy- 
geometry  are  presented  in  the  next  chapter. 


A.  Scatter-Free  Shield  Penetration 

This  test  problem  consists  of  a  uniform  ,  isotropic  flux  incident 
upon  the  left  of  a  shield,  with  vacuum  on  the  right  of  the  shield.  The 
shield  is  20  mean-f ree-paths  (mfp)  thick  and  is  purely  absorptive.  The 
exact  solutions  for  the  directional  flux  and  current  within  the  shield 
are  the  analytic  functions: 


F(y,x)  =  Fq  e 


-o  x  At 


Jx(U,x)  =  y Fq  e 


-  o  x  4i 


(IX-1) 


(IX-2) 


Integrating  equation  (IX-2)  over  jj  at  ax  =  20  gives  the  current 
penetrating  the  shield: 


Jx  =  0.5  E3(20)  =  4.50456  E-ll 


(IX-3) 


where  Fq  =  1  and  E^  is  an  exponential  integral. 

The  problem  was  solved  using  S„  with  equal  weight  quadrature  and 

N 

with  Gaussian  quadrature,  and  using  L  with  equal  weight  elements  and 

N 

Gaussian  weight  elements  for  various  angular  meshes  and  space  schemes. 
In  each  case,  40  space  cells  of  1/2  mfp  thickness  were  used.  For  this 
problem,  since  Q  =  0  in  every  cell,  the  step  characteristic  (SC)  and 
linear  characteristic  (LC)  methods  are  identical. 


VJll 


m 


Table  IX-1  presents  the  error  ratio  of  the  current  penetrating  the 
shield  for  the  various  schemes.  Error  ratio  is  a  computationally 
convenient  measure  of  error,  here  defined  as 


Error  Ratio  (x,x  >  =  lx  -  x  /  |x  +  x  | 

exact  1  ovart  1  ovarf 1 


exact 1 


exact 1 


(IJC-4) 


where  x  is  the  reference  or  benchmark  value.  This  has  the  advantage 
exact 

over  percentage  that  if  x  is  too  large  by  a  factor  of  2,  the  relative 
error  is  100%  but  if  x  is  too  small  by  a  factor  of  2,  it  is  -50%, 
whereas  the  error  ratio  is  1/3  in  each  case.  For  small  errors,  the  error 
ratio  is  approximately  half  the  conventional  relative  error.  An  error 
ratio  of  .05  corresponds  tc  approximately  10%,  for  example. 


Method  N:  2  4  6  6  12  24 

Three-Point-Gaussian  Discrete  Elements: 


Step/SC 

'ew 

.556 

.303 

.210 

.149 

Step/SC 

‘GW 

.556 

.220 

.0846 

.0407 

.0191 

.00504 

SC/SC 

EW 

.395 

.0488 

.00880 

SC/SC 

‘gw 

.395 

.0106 

.00122 

.000506 

Simpson's  Rule 

Discrete 

Elements: 

SC/SC 

‘ew 

.578 

.310 

.146 

.0670 

.0169 

sc/sc 

GW 

.578 

.161 

.0223 

.00930 

.00208 

.00012 

Discrete 

Ordinates : 

.  sc 

EW 

1.000 

.978 

.791 

.553 

.273 

.0696 

.  sc 
“g  q 

1.000 

.572 

.0759 

.00138 

.000007 

0 

N: 

8 

16 

24 

_32 

48 

SC 

‘ew 

.553 

.156 

.0696 

.0392 

.0174 

.  sc 
’gq 

.00138 

.000007 

.0000004 

0 

0 

Table  IX-1:  Error  Ratio  of  Shield  Penetration  Current 


1.  Step  vs .  SC  for  Auxiliary  Fluxes 

A  comparison  of  the  LsteP/sc  results  with  the  lSC/sc  results 
indicates  the  superiority  of  using  the  higher  order  method  for  the 
auxiliary  calculations.  Even  if  allowance  is  made  for  the  lower  cost  of 
the  step  method,  SC  is  more  effective.  For  example,  assuming  the  step 

method  costs  1/4  as  much  as  the  SC  method,  then  with  three  auxiliary 

Step/SC  SC/SC 

fluxes  per  element,  L„„t,  ^  costs  about  the  same  as  L_  ,  ,  but 

8EW  4EW 

has  three  times  as  much  error.  For  larger  N,  the  disparity  in 
performance  is  even  greater. 

2.  Simpson ' s  Rule  vs .  Gauss-Legendre  (3-point)  Quadrature 

SC/SC  sc/sc 

A  comparison  of  the  G3-L  and  SR-L  results  shows 

consistently  better  results  for  the  G3  method.  Even  if  the  SR  scheme 
reuses  the  overlapping  auxiliary  fluxes,  and  so  costs  only  3/4  as  much 
as  the  G3  scheme,  it  still  cannot  compete.  For  example,  the  G3-L6EW  and 
the  SR-L8ew  cost  about  the  same,  but  the  SR  scheme  has  over  7  times 
larger  error. 

3.  Equal  Weight  vs .  Gaussian  Weight  Discrete  Elements 

The  Gaussian  weight  discrete  elements  outperforms  the  equal 
weight  quadrature  mesh  for  this  problem.  The  reason  is  that  the  flux 
penetrating  completely  through  a  20  mfp  thickness  of  non-scattering 
absorber  is  strongly  forward-biased,  so  that  the  element (s)  closest  to 
!■!  =  1  carry  all  the  information.  With  Gaussian  weights,  the  discrete 
elements  are  crowded  toward  the  poles  and  sample  this  information  more 
effectively. 

4.  Discrete  Elements  vs .  Discrete  Ordinates 
When  compared  on  the  basis  of  number  of  elements/ordinates, 

SC/SC  schemes  have  smaller  error  than  the  SSC  schemes,  but  also 


the  G3-L 


have  slower  convergence  as  N  is  increased.  For  large  enough  N,  it  is 


expected  that  SN  will  be  the  more  accurate.  These  discrete  element 

methods,  however,  have  roughly  four  times  the  cost  per  element  of  the 

discrete  ordinates  schemes.  Comparing  S  to  L  ,  etc.,  the  Gauss-Legendre 

8  2 

quadrature  discrete  ordinates  method  is  quite  clearly  the  most  cost 
effective  scheme  tested,  for  this  problem.  The  bottom  rows  of  table  IX-1 
facilitate  this  comparison.  This  is  not  necessarily  a  representative 
comparison,  however,  because  the  problem  is  ideally  suited  for  the  SN 
quadrature  in  two  ways.  First,  the  angular  flux  distribution  which 
seeks  to  represent  with  a  high-order  polynomial  is  in  fact  an  analytic 
function  well  suited  to  such  a  representation,  namely,  the  function  in 
equation  (IX-1).  Secondly,  the  SC  spatial  quadrature  is  exact  (except 
for  rounding  error  accumulation)  for  this  source-free,  scatter-free 
medium,  so  that  the  data  points  used  to  obtain  the  high-order  fit  have 
essentially  no  error. 

B.  Two- Reg ion  Problem 

This  problem  consists  of  a  strongly  scattering  source  region  with 
weakly  scattering  shielding  and  vacuum  outer  boundaries.  Only  half  of 
the  system  is  solved,  using  a  symmetry  boundary  on  the  left  and  vacuum 
on  the  right.  The  problem  parameters  are  summarized  in  table  IX-2. 


Region 

Source  c 

Region  Width 

#  of  Cells 

Cell  Width 

(#/cm^/sec) 

(cm) 

Source 

1  0.5 

5 

40 

0.125 

Shield 

0  0.1 

15 

120 

0.125 

Table  IX-2:  Parameters  for  Two  Region  Test  Problem 


1 .  Error  "Norms 


Two  error  norms  are  used  in  each  region.  (The  term  "norm"  is 
used  loosely  hi  re;  no  attempt  is  made  to  demonstrate  that  these  are  true 
norms  in  the  function-analytic  sense.)  One  is  an  integral  norm;  the 
error  ratio  of  the  integral  of  the  scalar  flux  over  the  region.  This  is 
a  measure  of  the  ability  of  the  method  to  compute  reaction  rates  or 
related  eigenvalues.  The  other  is  a  pointwise  norm:  the  error  ratio  of 
the  scalar  flux  as  computed  for  each  spatial  cell,  averaged  over  the 
cells  of  the  region.  Tables  (IX-3)  through  (IX-6)  present  these  results. 


Method  Quadrature :  2 

_4 

6 

8 

48 

Three-Point 

-Gaussian 

Discrete  Elements: 

Step/DDF 

lew 

.000112 

.000058 

.000044 

.000038 

Step/LC 

lew 

.000097 

.000040 

.000025 

.000019 

Step/LC 

lgw 

.000097 

.000054 

.000035 

.000026 

LC/LC 

lew 

.000021 

.000012 

.000011 

.000010 

LC/LC 

lgw 

.000021 

.000014 

.000012 

.000011 

Simpson' 

s  Rule  Discrete  Elements: 

LC/LC 

lew 

.000103 

.000004 

.000006 

.000009 

Gaussian  Quadrature  Discrete  Ordinates: 

sstep 

.00619 

.00200 

.000983 

.000639 

sDDF 

.00534 

.00125 

.000286 

.000025 

ssc 

.00537 

.00128 

.000329 

.000022 

sLC 

.00535 

.00126 

.000553 

.000307 

Benchmark 

Quadrature:  8 

_16 

24 

sLC 

.000307 

.000071 

.000027 

Table  IX-3:  Error  Ratio  of  Scalar  Flux  Integrated  over  Source  Region 


2.  Step/LC  vs.  Step/DDF 

,ie  least  expensive  discrete  element  scheme,  per  cell  per 

Step/DDF 

element,  which  might  produce  reasonable  accuracy,  is  the  L 
scheme.  The  step  auxiliary  calculations  are  smooth  and  the  DDF  main 
calculations  are  relatively  accurate.  Comparison  of  the  results  with 
those  of  the  L  F  scheme  shows  the  latter  to  be  more  accurate  on  an 
equal  N  basis  and  to  have  faster  convergence  in  N,  but  (allowing  for 
twice  the  cost  per  element)  the  two  are  about  equally  cost-effective. 


Table  IX-4:  Error  Ratio  of  Scalar  Flux  Integrated  over  Shield  Region 


3.  Equal  Weight  vs.  Gaussian  Weight  Elements 


For  both  integrated  and  pointwise  errors  and  in  both  source 
and  shield  regions,  LgW  is  more  accurate,  and  hence  more  cost-effective, 
than  L  .  The  only  exception  in  the  four  tables  is  the  L^  average  error 
ratio  in  the  shield  region.  This  observation  holds  for  both  the  Step/LC 
and  the  LC/LC  spatial  quadratures.  The  introduction  of  scatter  and  an 
interior  material  boundary  make  this  problem  more  realistic  than  the 
non-scattering  shield,  and  show  the  validity  of  equal  weight  composite 
quadrature  for  less  well-behaved  fluxes. 


Method  Quadrature :  2 

4 

6 

8 

48 

Three-Point -Gaussian 

Discrete  Elements: 

Step/DDF 

EW 

.00106 

.000540 

.000392 

.000356 

Step/LC 

EW 

.000816 

.000276 

.000142 

.000087 

Step/LC 

GW 

.000816 

.000441 

.000253 

.000162 

_  LC/LC 

lew 

.000314 

.000073 

.000057 

.000046 

LC/LC 

GW 

.000314 

.000109 

.000068 

.000061 

Simpson ' s 

Rule  Discrete  Elements: 

LC/LC 

lew 

.00169 

.000433 

.000186 

.000^03 

Gaussian  Quad 

~ature  1 

Discrete  Ordinates: 

sstep 

.00761 

.00258 

.00211 

.00259 

sDDF 

.00855 

.00241 

.000780 

.000303 

ssc 

.00850 

.00232 

.000639 

.000056 

sLC 

.00846 

.00227 

.00104 

.000597 

Benchmark 

Quadrature:  8 

16 

24 

sLC 

.000597 

.000149 

.000058 

4.  LC/LC  vs.  Step/LC 


For  both  integrated  and  pointwise  errors,  and  in  both  source 

and  shield  regions,  and  even  allowing  for  a  2:1  or  greater  cost  ratio, 

LC/LC  .  Step/LC 

L  is  more  accurate  and  more  cost-effective  than  L  .  Similar 

SC /SC 

comparisons  with  L  (data  not  shown)  reveal  the  same  trend:  the 

higher-order  the  spatial  quadrature,  the  more  accurate  and  cost- 
effective  the  discrete  element  method.  The  use  of  LC  for  both  auxiliary 
and  main  fluxes  can  be  even  more  efficient  if  its  ability  to  use  a 
coarser  spatial  mesh  is  used  to  advantage. 


Method  Quadrature :  2 

6 

8 

48 

Three-Point 

-Gaussian 

Discrete  Elements: 

Step/DDF 

lew 

.0474 

.0270 

.0188 

.0139 

Step/LC 

lew 

.0404 

.0203 

.0124 

.00770 

Step/LC 

GW 

.0404 

.0130 

.0446 

.00287 

T  LC/LC 

lew 

.0340 

.000914 

.000094 

.000027 

_  LC/LC 

lgw 

.0340 

.000483 

.000123 

.000033 

Simpson' 

s  Rule  Discrete  Elements: 

_  LC/LC 
lew 

.203 

.0425 

.0107 

.00360 

Gaussian  Quadrature  Discrete  Ordinates: 

s5tep 

.601 

.185 

.223 

.224 

S°DF 

.707 

.0785 

.00700 

.00593 

ssc 

.704 

.0737 

.000915 

.000260 

s1* 

.704 

.0738 

.00324 

.00113 

Benchmark 

Quadrature :  8 

16 

24 

s^ 

.00113 

.000231 

.000088 

Table  IX-6:  Average  Error  Ratio  of  Scalar  Flux  in  the  Shield  Region 
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5.  £2  vs .  SR  Element  Quadrature 

For  pointwise  errors,  in  both  source  and  shield  regions,  the 

three-point  Gauss-Legendre  quadrature  provides  clearly  superior 

performance.  Even  allowing  for  a  2:1  cost  ratio  (the  actual  ratio  is 

LC/LC 

closer  to  1.2:1),  the  G3-L  scheme  is  more  cost-effective  than  the 

EW 

LC/LC 

SR-L„,,  scheme. 

EW 

For  the  integrated  flux  errors,  the  evidence  is  not  so  clear-cut. 

In  the  source  region,  for  example,  both  methods  converge  toward  an  error 

ratio  of  0.000010,  indicating  this  much  residual  error  in  the  benchmark 
Lc 

(S48  )  solution.  But  they  converge  from  opposite  directions.  The  G3-L^ 

is  nearly  converged,  at  0.000012,  and  is  actually  closer  than  the  SR-L-. 

o 

which  is  also  nearly  converged,  at  0.000004.  However,  neither  of  these 
errors  is  large  enough  for  meaningful  comparisons.  Similar  logic  applies 
to  the  data  for  the  shield  region. 

6.  Discrete  Elements  vs.  Discrete  Ordinates 

Comparison  of  the  best  discrete  elements  method  tested,  namely 
LC/LC 

G3-L  ,  with  the  best  discrete  ordinates  quadrature  tested,  S  with 

EW  N 

Gauss-Legendre  quadrature,  shows  that  the  L  method  is  very  much  more 

N 

accurate,  for  the  same  N.  On  a  cost  basis,  allowing  for  a  4:1  cost 

ratio,  the  bottom  row  of  each  table  can  be  compared  directly  with  the  L_, 

rows.  In  oerms  of  the  integrated  flux  errors,  the  L„  method  is  the 

better  performer.  Even  L.  has  less  error  than  S,„,  and  at  1/4  the 
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computational  cost.  In  terms  of  the  average  errors  in  the  source  region, 
L^  is  somewhat  more  cost-effective.  Only  in  terms  of  the  average  errors 
in  the  shield  region  is  SN  the  more  accurate  for  a  given  cost.  This 
occurs  for  reasons  similar  to  those  which  applied  to  the  non-scattering, 
non-source  shield  penetration  problem. 
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It  might  be  argued  that  a  double-range  Gauss-Legendre  quadrature 
would  improve  the  performance  of  the  discrete  ordinates  method  for  this 
problem  by  modeling  the  flux  discontinuity  at  u  =  0  exactly,  as  with 
the  double  PN  method.  With  this  quadrature,  SN  could  possibly  be  more 
cost-effective  than  L„.  However,  Carlson  and  Lathrop  have  reported  that 
the  single-range  quadrature  is  more  accurate  for  problems  with  optically 
thick  regions,  while  double-range  is  better  for  problems  with  many  thin 
regions  [Ref.  4:34].  The  regions  in  this  problem  are  optically  thick  (10 
and  15  mfp),  so  that  the  comparison  with  single-range  SN  should  be  a 
fair  one. 

C.  Conclusions 

This  chapter  has  presented  the  results  of  numerical  testing  of 
various  discrete  element  and  discrete  ordinate  schemes  in  slab  geometry. 
As  a  proof  of  concept,  two  test  cases  were  used.  While  more  extensive 
testing  would  be  needed  to  support  categorical  conclusions,  the  data 
presented  do  indicate  the  following  trends. 

1  -  Equal  weight  discrete  elements  are  generally  more  accurate  than 
Gauss-Legendre  weight  elements,  for  realistic  problems. 

2  -  Three-point  Gauss-Legendre  element  quadrature  is  superior  to 
Simpson's  rule. 

3  -  L  method  can  be  more  cost-effective  than  Gaussian  S  when  the 

N  N 

linear  characteristic  space  scheme  is  used  for  both  the  auxiliary  and 
main  fluxes,  and  the  three-point  Gauss-Legendre  rule  is  used  for  the 
element  quadrature. 

4  -  The  expectation  that  "the  harder  the  problem,  the  more 
advantageous  the  L  method"  is  supported  by  this  evidence. 

N 


5  -  Convergence  and  accuracy  of  are  degraded  by  use  of 
inaccurate  or,  especially,  non-smooth  space  schemes  for  the  auxiliary 
fluxes.  Ln  with  step  characteristic  or  linear  characteristic  auxiliary 
quadrature  converged  in  the  same  number  of  iterations  as  S^,  however. 

The  above  tests  in  one  dimension  implicitly  avoided  consideration 
of  ray  effects,  since  there  are  no  ray  effects,  per  se,  in  one 
dimensional  (time  independent)  problems.  The  next  chapter  uses  two  test 
problems,  one  of  which  is  similar  to  the  ray  effects  problem  of  chapter 
V.  An  objective  of  those  tests  is  to  determine  the  applicability  in  two 


dimensions  of  the  observations  made  here  in  one  dimension. 


M 
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This  chapter  compares  the  performance  of  and  schemes  for  two 
problems  in  two-dimensional  Cartesian  geometry.  The  first  problem  is  a 
variation  of  the  square  source  in  a  square  absorber  problem  introduced 
in  chapter  VI.  A  review  of  results  from  the  literature  provides 


qualitative  comparisons  of  performance  in  ameliorating  ray  effects.  The 
second  problem  models  a  shielded  plane  source  with  a  vacuum  duct  in  the 
shield,  and  is  used  for  quantitative  as  well  as  qualitative  comparisons 
of  performance  on  a  difficult  problem. 


A.  Ray  Effect  Problem  and  Comparison  of  Results 

A  test  problem  which  has  been  used  in  the  literature  for  evaluation 


of  ray  effect  is  inset  in  figure  X-l.  The  problem  is  similar  to  the  one 
used  in  chapter  VI,  but  with  the  source  square  increased  in  size  and 
with  scatter  included.  The  results  used  for  comparison  are  the  cell- 
average  scalar  fluxes  along  the  top  row  (or  up  the  right  column,  by 
symmetry)  of  a  30  cell  by  30  cell  spatial  mesh.  This  mesh  is  nearly 
twice  as  fine  as  that  used  in  chapter  VI,  so  the  quasi-ray  effect  should 


be  reduced. 


1.  Conventional  Discrete  Ordinates 


Lathrop  used  this  problem  with  several  angular  quadrature 


schemes.  Figures  X-l  through  X-4  are  taken  directly  from  reference  8  for 

qualitative  comparisons.  Figure  X-l  shows  the  results  for  S2,  S^,  and 

S, , ,  using  conventional  (totally-symmetric)  quadrature  using  TWOTRAN 
lo 

(and  presumably,  diamond  difference  spatial  quadrature).  The  substantial 
wobbles  in  the  curves  are  described  by  Lathrop  as  ray  effect. 


m 


2.  Rotationally  Synunetric  Discrete  Ordinates 

Figure  X-2  shows  the  results  using  discrete  ordinates  with  a 
product  quadrature  set  consisting  of  a  single  row  of  directions  at  a 
fixed  latitude  (polar  angle),  equally  spaced  in  azimuth  and  equally 
weighted,  i.e.,  the  "rotationally  invariant"  quadratures.  In  Lathrop's 
notation,  these  are  R  ,  R  ,  R  ,  and  R  .  These  show  decreased,  but 
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still  significant,  ray  effect. 

3.  "Consistent"  Discrete  Ordinates 

Figure  X-3  compares  the  S_.  result  with  that  of  Lathrop's  CS_ 

io  2 

method.  This  "consistent  S^"  method  [Ref.  10}  used  explicit  coupling  of 
the  directional  flux  distribution  to  the  spatial  quadrature  (as  does 
discrete  elements)  derived  by  a  first-order  Taylor  expansion  of  the  flux 
about  the  SN  quadrature  angles.  This  scheme  shows  a  vast  improvement 
over  the  S2  scheme  (to  which  it  degenerates  if  the  coupling  is  omitted), 
and  might  seem  quite  promising.  Its  limitation  seems  to  be  that  the 
"consistent"  treatment  is  an  approach  rather  than  a  numerical  method,  in 
that  to  increase  either  the  set  of  quadrature  angles  or  the  order  of  the 
coupling  terms,  or  both,  requires  a  complete  rederivation  of  the 
equations  to  be  solved.  With  the  discrete  elements  method,  however, 
increasing  the  quadrature  set  is  a  strictly  mechanical  procedure,  and 
increasing  the  coupling  order  requires  only  the  use  of  a  different 
element  quadrature  rule. 

4.  Spherical-Harmonics-Like  Discrete  Ordinates 

In  any  case,  Lathrop  abandoned  the  CS2  scheme  in  favor  of 
"spherical-harmonics-like"  discrete  ordinates.  This  method  uses 


fictitious  sources  to  cause  the  discrete  ordinates  equations  to  conserve 


the  same  moments  as  the  spherical  harmonics  method  [Ref.  8).  Results 


with  these  schemes,  S  -•»  P, ,  S.-->  P. ,  S  -■»  P,,  S_— >  P_,  are  shown  in 

2  14  Id  38  D 

figure  X-4.  They  show  no  "ray  effect"  in  that  they  are  shaped  correctly, 
but  they  have  angular  quadrature  truncation  error  which  systematically 
affects  the  level  rather  than  the  shape.  The  convergence  in  amplitude 
(vs.  N)  is  relatively  slow,  as  with  the  PN  method.  The  method  also  has 
substantially  increased  computational  cost. 

5.  Compatible  Product  Quadrature  Discrete  Ordinates 

Abu-Shumays  [Ref.  2]  investigated  the  possibilities  of 
improved  performance  through  use  of  specialized  "compatible"  product 
quadrature  sets  (with  diamond  difference  spatial  quadrature).  He 
concluded  that  the  Gauss-Christof fel  polar  /  quadruple-range  azimuthal 
product  quadrature  has  more  accuracy  and  less  ray  effect,  for  this 
problem,  than  any  of  the  other  quadratures  he  investigated.  The  results 
for  this  quadrature,  S,_  and  S  .  are  shown  in  Figure  X-5,  taken 

directly  from  reference  2. 

6.  Discrete  Elements 

A  variety  of  discrete  element  schemes  were  tested  for  this 
problem  using  a  16  by  16  spatial  grid.  For  the  auxiliary  fluxes,  the 
step  method  was  found  to  be  unacceptable  in  accuracy,  while  the  diamond 
difference  method  (with  or  without  negative  fixups)  was  inaccurate  and 
slow  to  converge  (L„_  .),  or  prevented  convergence  (L  ).  with  the  step 
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characteristic  method  for  the  auxiliary  fluxes  and  either  diamond 
difference  or  step  characteristic  for  the  main  fluxes,  convergence  was 
as  fast  as  for  the  SN  method.  As  in  the  one-dimensional  test  problems, 
the  Gauss-Legendre  three-point  element  quadrature  was  generally  more 


accurate  than  the  Simpson's  Rule.  The  Newton-Cotes  rule  was  about  as 


accurate  as  the  Gauss-Legendre,  but  at  greater  execution  cost. 


The  results  of  two  calculations  with  the  30  by  30  spatial  grid  are 

SC /DDF 

shown  in  figures  X-6  and  X-7.  The  schemes  used  were  G3-L  and 
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G3-L  .  The  SC/DDF  solution  shows  a  small  amplitude  wiggle  which 
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can  be  attributed  to  the  spatial  quadrature,  since  it  is  not  present  in 
the  SC/SC  solution,  and  so  is  quasi-ray  effect.  The  SC/SC  solution  also 
shows  a  quasi-ray  effect,  in  that  there  is  a  slight  break  in  slope  at 
x  *  1  which  corresponds  to  the  edge  of  the  source  region  (as  was  seen 
in  chapter  VI).  These  quasi-ray  effects  were  reduced  by  refining  the 
spatial  grid  from  16  by  16  to  32  by  32.  The  SC/SC  curve  falls  between 
the  S  ->  P,  and  S  ->  Pc  curves  of  figure  X-4. 

Comparing  figure  X-6  with  figures  X-l  through  X-5,  it  may  be 
concluded  that,  for  this  problem,  the  discrete  element  method,  with  as 
few  as  three  elements  in  azimuth  per  azimuthal  quadrant,  ameliorates  ray 


effect  as  well  as  the  best  of  previous  methods. 


Fig.  X-l:  Ray  Effects:  Conventional  SN  [Ref.  8:257) 


Rtflectivt 

Boundary 


j  +  i  f- 


Vacuum  Duct  Problem:  Qualitative  Comparison  of  Results 


The  second  test  case  is  a  duct  problem  which  is  difficult  for 
discrete  ordinates.  A  quantitative  comparison  of  L  and  S  is  made. 
Streaming  ducts  in  shields  are  a  common  design  problem  for  which 
empirical  thumb-rules  are  often  used.  Where  better  solutions  are  needed, 
Monte  Carlo  methods  provide  accurate  but  expensive  answers.  This  section 
presents  graphs  of  numerical  results  and  qualitative  comparison;  the 
next  section  presents  quantitative  comparisons. 

1.  The  Test  Problem 

Figure  X-8  defines  a  test  problem  consisting  of  a  thin  source 
region  along  the  bottom  of  a  shield  with  vacuum  boundaries  along  the 
top,  right,  and  bottom.  The  left  side  of  the  problem  is  a  symmetry 
boundary,  with  a  deep,  narrow  vacuum  duct,  or  channel,  along  the  edge. 
It  is  expected  that  the  major  leakaqe  through  the  shield  will  be  by 
streaming  up  the  duct.  This  particular  problem  was  selected  as  an 
idealized  representation  of  an  access  port  in  a  fusion  reactor  design. 
Such  access  ports  are  required  for  plasma  injection,  charged  particle 


heating  beams,  laser  instrumentation,  etc. 


2.  Monte  Carlo  Benchmark  Solution 

In  order  to  accurately  assess  the  performance  of  both  the 
discrete  elements  and  discrete  ordinates  methods,  a  benchmark  solution 


known  to  be  both  accurate  and  free  of  ray  effect  was  required.  The  Monte 

Carlo  solution  used  for  this  purpose  is  shown  in  figures  X-9  through 
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X-13.  The  solution  entailed  the  simulation  of  10  particles.  Data 

recorded  included  total  leakage,  total  absorption,  and  leakage  through 

each  of  16  equal  segments  along  the  top  edge  of  the  problem, 

corresponding  to  the  16  cell  by  16  cell  spatial  grid  used  for  the 

discrete  elements  and  discrete  ordinates  calculations.  The  relative 

variance  of  the  results  was  approximately  1%  at  the  one-sigma  level  of 

confidence  for  the  least  accurate  of  the  16  top  face  leakages.  This 

solution  was  provided  by  W.  T.  Urban  of  Los  Alamos  National  Laboratory 

in  support  of  this  research. 

3.  Conventional  Discrete  Ordinates 

Conventional  S„  quadrature  solutions  using  various  spatial 
N 

quadrature  schemes  were  provided  by  E.  W.  Larsen  in  support  of  this 

research.  Diamond  difference  solutions  are  shown  in  figure  X-9.  There 

DDF 

are  two  striking  observations  to  be  made  about  SN  :  the  results  are 

not  good,  and  they  don’t  improve  much  as  N  is  increased.  The  currents 

are  in  error  by  factors  of  two  to  four  over  most  of  the  top  face  of  the 

problem,  even  for  S..,  which  uses  36  quadrature  directions  per  octant. 

16 

The  step  characteristic  solutions  are  smoother  and  more  accurate  in 
level,  and  are  shown  in  figure  X-10.  Larsen  and  Alcouffe  have  developed 
a  linear  characteristic  spatial  quadrature  method  for  xy-geometry 
[Ref.  5J.  Results  with  this  quadrature,  shown  in  figure  X-ll,  have  some 


improvement  over  the  step  characteristic  method. 
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4.  Product  Quadrature  Discrete  Ordinates 


Figure  X-12  shows  the  results  of  a  discrete  ordinates 

calculation  using  a  Gauss-Christof fel  /  rotationally  symmetric  product 

SC 

quadrature  and  step  characteristic  spatial  quadrature:  S^_  -  .  This  is 
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the  angular  quadrature  underlying  the  discrete  elements  results 
considered  next.  These  results  are  greatly  improved  over  the 
conventional  S„  results,  both  in  overall  level  and  in  shape.  The 
conventional  SN  gave  a  ratio  of  less  than  4:1  for  the  leakage  in  the 
duct  to  that  on  the  right  half  of  the  top  of  the  shield,  even  for  S.  . 
The  product  quadrature  gives  a  ratio  of  about  6:1,  while  the  benchmark 
ratio  is  about  12:1.  In  the  next  section,  a  "level  norm"  and  a  "shape 
norm"  are  introduced  to  properly  quantify  these  observations. 

5.  Discrete  Elements 

Figure  X-13  shows  the  results  for  the  discrete  elements 
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method:  G3-L__  ,  ,  with  16  by  16  and  32  by  32  spatial  meshes.  The 

ratio  of  duct  current  to  that  on  the  right  is  about  11:1,  much  improved 
over  all  the  discrete  ordinates  methods.  Also,  the  leakage  is  monotone 
decreasing  from  left  to  right,  as  expected  from  the  problem  symmetry. 
None  of  the  discrete  ordinates  methods  were  correct  in  this  regard. 
Deficiencies  of  this  discrete  element  solution  are: 

a  -  The  leakage  through  the  duct  is  overestimated  by  20%-30% 
b  -  There  is  a  slight  hump  to  the  curve  at  about  x  ■  1.5  cm 
c  -  The  leakage  through  the  shield  away  from  the  duct  is 
overestimated  by  about  15% 

Item  b  is  a  ray  effect  and  is  eliminated  in  the  4  solution.  The 
other  two  deficiencies  are  at  least  partially  quasi-ray  effect  and  are 
reduced  by  the  refined  (32  by  32)  spatial  mesh. 
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C.  Vacuum  Duct  Problem:  Quantitative  Comparison  of  Results 


This  section  makes  a  quantitative  comparison  of  the  accuracy  of 
discrete  elements  solutions  with  that  of  discrete  ordinates  solutions, 
for  the  problem  presented  in  the  previous  section.  As  in  chapter  IX, 


informal  "norms"  are  used  for  this  purpose.  Two  such  measures  of  error 
are  used.  One  is  sensitive  to  the  over-all  amplitude  of  the  top  face 
leakage.  This  "level  norm"  is  the  error  ratio  (defined  by  equation  IX-4) 
of  the  total  leakage  through  the  top  face  of  the  problem,  i.e.,  through 


the  duct  and  through  the  shield,  with  respect  to  the  Monte  Carlo 
benchmark  solution.  The  table  indicates  the  sign  of  the  error:  positive 
for  overestimated  leakage,  negative  for  underestimated  leakage.  The 
other  measure  of  error  is  sensitive  to  the  relative  distribution  of 
leakage  through  the  top  face.  This  "shape  norm"  is  the  average  of  the 
error  ratios  of  the  16  normalized  cell  currents.  The  currents  are 
normalized  to  give  the  same  total  leakage  as  the  benchmark  solution,  so 
that  the  norm  will  be  sensitive  to  shape  but  not  level.  The  performance 
of  various  schemes  tested  is  presented  in  table  X-l,  in  terms  of  these 
norms.  Data  is  presented  for  the  step  characteristic  spatial  quadrature, 
since  the  SC/DDF  results  show  poor  performance  in  the  shape  norm  as  a 


result  of  spatial  oscillations  (quasi-ray  effect).  S  results  are 
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also  included.  Observations  based  on  this  data  are  presented  next. 

1.  Convergence  of  Conventional  SN 

Refining  the  angular  quadrature  from  S_  to  S,c  provided 

o  lb 

negligible  quantitative  improvement  in  either  shape  or  level.  Refining 
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the  spatial  quadrature  from  to  ,  however,  improved  the  level 

but  not  the  shape.  The  conventional  S  angular  quadrature  is  essentially 
converged,  and  with  poor  results. 
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2.  Convergence  in  Polar  Quadrature 

Increasing  the  accuracy  of  polar  quadrature  by  increasing  K 

for  product  quadratures  should  improve  level  but  have  little  effect  on 

shape.  This  effect  is  seen  in  a  comparison  of  S__  and  S  ,  both 
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of  which  are  essentially  converged  in  azimuth.  Also,  the  hybrid  discrete 

element  method,  G3-L  is  significantly  more  accurate  in  level  than 
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the  corresponding  G3-L  method,  but  only  a  little  better  in  shape, 
indicating  its  superiority  in  polar  quadrature.  The  hybrid  method  is 
more  accurate  and  less  expensive. 

3.  Convergence  in  Azimuthal  Quadrature 

Each  of  the  product  quadratures  S0_  and  S,_  ,  use  36 
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ordinates  per  octant  and  thus  have  the  same  computational  cost  as  S^. 
Because  they  apply  most  of  their  effort  to  the  difficult  azimuthal 
quadrature,  and  because  they  use  a  composite  quadrature  scheme 
(composite  midpoint  rule)  to  treat  the  ill-behaved  azimuthal  flux 
variations,  they  greatly  outperform  the  conventional  quadrature.  The 
azimuthal  quadrature  is  converged  at  a  shape  norm  of  about  .035-.038  for 
the  various  product  quadrature  discrete  elements  and  discrete  ordinates 
schemes  tested,  for  a  16  cell  by  16  cell  spatial  grid.  Thus,  to  the 
accuracy  of  this  spatial  quadrature,  the  G3-L0_  ,  scheme  is  essentially 
converged  in  azimuthal  quadrature.  The  Simpson's  rule  method,  however, 
is  less  effective,  not  being  converged  in  shape  at  SR-L__  The  Newton- 
Cotes  rule  method  is  as  accurate  as  the  Gauss-Legendre  method,  but  costs 


more  to  compute. 


The  composite  Gauss-Christof fel  /  quadruple-range  quadrature, 

S  was  evaluated  for  both  space  mesh  sizes.  The  level  error 

decreased  with  the  refined  spatial  grid,  but  the  shape  error  increased. 

This  inconsistent  behavior,  together  with  the  average  19%  pointwise 

error,  indicates  that  the  nearly  exact  total  leakage  for  the  method  is 

coincidental,  as  for  the  .  . 
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The  improvement  in  the  G3-L__  method  resulting  from  the  spatial 
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grid  refinement  was  more  significant  than  that  resulting  from  a 

refinement  of  the  angular  grid  to  G3-L  (still  on  a  16x16  spatial 
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grid).  This  improvement  was  both  in  level  and  in  shape.  The  reasons  for 
these  improvements  are  considered  below.  With  this  refined  spatial  grid, 
further  improvement  was  noted  by  then  refining  the  angular  mesh  as  well, 
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Since  the  errors  of  shape  and  level  remaining  in  the  discrete 
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elements  solutions,  such  as  G3-L  ,  are  dominated  by  the  spatial 

§  <3 

quadrature,  and  since  it  was  not  possible  to  refine  the  spatial  mesh  to 
convergence  within  the  scope  of  this  research,  precise  comparisons  of 
cost-effectiveness  of  the  various  schemes  could  not  be  made.  A  more 


important  consideration,  perhaps,  is  the  confidence  level  (in  an 
engineering  sense)  of  the  various  methods,  which  is  considered  next. 


5.  Confidence  of  Results 


The  foregoing  discussions  of  accuracy  and  convergence  have 
implicitly  presumed  that  the  benchmark  solution  is  known.  From  the 
viewpoint  of  practical  engineering,  however,  the  objective  of  numerical 
calculations  is  to  estimate  an  (unknown)  answer  and  its  confidence 
limits. 

Taking  this  view,  and  supposing  the  duct  problem  is  to  be  solved 
only  by  the  discrete  ordinates  method,  the  results  obtained  above  lead 
to  a  low  level  of  confidence.  Different  families  of  quadrature  sets 
apparently  converge  to  different  results,  both  in  shape  and  level.  The 
convergence  is  inconsistent  with  respect  to  the  two.  In  terms  of  the 
error  valleys  that  Lathrop  described,  the  dilemma  is  simply  stated:  "how 
is  one  to  know  which  solutions  are  in  the  valley?"  As  seen  with  the 
total  top-face  leakage  and  the  S„_  __  quadratures,  the  hills  on  one  side 
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may  be  positive  error  while  the  hills  on  the  other  side  of  the  valley 
are  negative  error.  The  properly  conservative  engineering  answer  is  that 
the  shield  leakage  is  the  median  of  the  collection  of  discrete  ordinates 
results,  to  within  a  factor  of  about  two  (for  this  example).  This  is  the 
reason  nuclear  power  reactor  shielding  is  typically  over-designed. 

If  this  same  engineering  problem  were  to  be  solved  by  discrete 

elements,  however,  the  consistent  convergence  in  polar,  azimuthal,  and 

spatial  quadratures,  toward  the  same  answer,  using  different  element 

quadrature  rules  and  different  quadrature  types  (L.  vs.  L  T )  allow  a 
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much  more  useful  conclusion  to  be  drawn.  The  result  can  be  estimated  to 
a  confidence  of  about  plus  or  minus  10%.  There  may  still  be  an  "error 
valley",  but  all  the  results  in  table  X-l  are  on  one  side  jf  it,  and  the 
hills  are  lower. 
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D.  Interaction  of  Spatial  and  Angular  Quadrature 

The  discrete  elements  solutions  for  the  vacuum  duct  problem  were 

not  fully  converged  in  spatial  quadrature.  As  a  consequence,  the 

auxiliary  fluxes  were  not  as  accurate  as  possible,  so  that  the  coupling 

of  angular  and  spatial  quadratures  is  not  as  effective  as  possible.  The 
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errors  of  the  G3-L  solution  were  observed  to  be  of  three  types: 
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1  -  The  leakage  through  the  duct  is  overestimated  by  20%-30% 

2  -  There  is  a  slight  hump  to  the  curve  at  about  x  =  1.5  cm 

3  -  The  leakage  through  the  shield  away  from  the  duct  is 
overestimated  by  about  15%. 

The  second  item  is  a  ray  effect  and  is  eliminated  by  refining  the 
angular  mesh,  but  the  other  two  are  quasi-ray  effects,  and  are  reduced 
by  refining  the  spatial  mesh.  The  objective  cc  this  section  is  to 
evaluate  the  causes  of  these  quasi-ray  effects  and  the  possibilty  of 
reducing  or  eliminating  them  through  the  use  of  a  higher  order  spatial 
quadrature . 

1 .  Leakage  Through  the  Shield 

Since  it  is  a  simpler  problem,  the  overestimation  of  the 
leakage  through  the  shield,  away  from  the  vacuum  duct,  is  considered 
first.  The  step  characteristic  scheme  assumes  a  flat  source  distribution 
within  each  spatial  cell.  This  source  term  is  Q  and  includes  scatter.  As 
the  flux  penetrates  a  shield,  it  should  be  monotone  decreasing,  so  that 
the  scatter-source  should  be  decreasing  across  each  cell  (in  the 
direction  of  shield  penetration).  But,  the  step  characteristic  scheme 
computes  the  cell-average  value,  redistributing  Q  uniformly  across  the 
cell.  This  is  a  systematic  error  which  consistently  transports  scattered 
particles  deeper  into  the  shield,  as  shown  in  figure  X-14A. 


A:  Step  Characteristic  B:  Linear  Characteristic 
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Fig.  X-14:  SC  and  LC  Shield  Penetration  Errors 

The  linear  characteristic  method  models  the  scatter-source  term,  Q, 
as  linearly  varying  across  each  spatial  cell  so  as  to  approximate  both 
its  average  and  its  first  (spatial)  moment.  This  eliminates,  to  the 
accuracy  of  the  approximation,  the  systematic  artificial  transport  of 
particles  into  the  shield,  as  shown  in  figure  X-14B.  It  is  expected, 
then,  that  use  of  the  LC  spatial  quadrature  with  the  discrete  elements 
method  would  substantially  reduce  the  overestimation  of  the  leakage 
through  the  shield,  away  from  the  streaming  duct. 

2.  Leakage  Through  the  Duct 

The  presence  of  the  duct  enhances  overall  leakage  through  the 
shield  in  two  ways.  Particles  that  originate  below  the  duct  (or  scatter 
below  the  duct)  may  stream  freely  through  the  shield,  provided  only  that 
they  start  in  the  right  direction.  This  is  illustrated  as  path  A  of 
figure  X-15.  Discrete  ordinates  is  inaccurate  for  this  problem  because 
few,  if  any,  of  the  quadrature  directions  represent  this  path.  Discrete 
elements,  however,  can  steer  the  streaming  into  this  path.  The  second 
effect  of  the  duct  is  to  decrease  the  optical  thickness  of  the  shield 
for  paths  which  cross  through  the  duct  and  ultimately  leak  out  through 
the  main  body  of  the  shield,  as  illustrated  by  path  B. 


Fig.  X-15:  Paths  Through  the  Vacuum  Duct 


The  discrete  elements  method  errs  in  a  different  way  in  handling 
path  B.  The  flux  in  an  angular  element  at  point  C  of  figure  X-15  would 
consist  of  mostly  path  A  type  particles  with  an  admixture  of  some  path  B 
type  particles,  and  would  be  represented  by  a  flux-weighted  average 
streaming  direction  such  as  D.  Assuming  analytic  spatial  quadrature  were 
somehow  employed,  the  particles  in  the  element  would  be  moved  precisely 
in  direction  D;  the  result  would  be  to  walk  the  B’s  across  the  duct  and 
the  A's  up  the  duct.  But,  with  step  characteristic  spatial  quadrature, 
this  is  not  the  result.  Since  direction  D  is  close  to  the  duct  axis  (y), 
and  hence  far  from  the  cell  diagonal,  the  method  is  numerically 
inaccurate  and  (as  shown  in  section  VI-C)  walks  the  flux,  like  a  chess 
rook,  straight  up  the  duct.  The  renormalizing  of  the  spatial 
distribution  along  the  cell  edge  to  a  flat  approximation  moves  the 


particles  (falsely)  to  the  left  at  each  cell  interface.  The  path  B 
particles,  in  effect,  become  trapped  in  the  duct  and  escape  absorption. 


The  linear  characteristic  scheme  would  mimimize  this  form  of  error 
by  representing  the  flux  along  each  cell  interface  with  a  linear  fit 
which,  at  least  approximately,  matches  the  average  and  the  first 
(spatial)  moment  along  the  cell  boundary.  This  eliminates  the  flattening 
process  which  moved  particles  leftward,  and  would  allow  the  B*s  to  cross 
the  duct  and  the  A's  to  stream  up  the  duct.  The  B's  would  not  arrive  at 
the  top  of  the  duct,  so  that  the  overestimation  of  the  duct  leakage  by 
the  step  characteristic  scheme  should  be  largely  eliminated.  Since  the 
B's  would  be  exposed  to  further  absorption  upon  reentering  the  shield 
material,  only  some  would  survive  to  leak  through  the  shield,  so  that 
the  total  leakage  would  be  reduced.  This  same  feature  would  reduce  or 
eliminate  the  quasi-ray  effect  of  the  SC  scheme  seen  in  the  square-in-a- 
square  problems  in  sections  VI-C  and  X-A. 

The  linear  characteristic  scheme  is  not  merely  higher  order  than 
the  step  characteristic  scheme,  but  it  also  corrects  the  specific 
deficiencies  which  are  the  source  of  the  errors  observed  in  the  discrete 
elements  solutions.  Thus,  the  linear  characteristic  method  is  expected 
to  correct  both  the  level  and  shape  errors  of  the  step  characteristic 

method,  when  used  for  both  the  main  and  auxiliary  flux  calculations  of 
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the  discrete  element  method.  A  demonstration  of  the  validity  of  this 
logic  by  actually  programming  the  combined  linear  characteristic  / 
discrete  elements  method  was  beyond  the  scope  of  this  project,  but 
refining  the  spatial  mesh  ameliorates  the  deficiencies  of  the  SC  scheme, 
and  was  observed  to  improve  both  the  shape  and  level  of  the  results  in 
the  manner  anticipated  from  this  analysis. 

*  The  linear  nodal  method,  recently  developed  by  Walters  and  O'Dell 
(Ref.  13]  should  also  be  appropriate  for  this  application. 


E.  Conclusions 


The  conclusions  derived  from  the  one-dimensional  test  cases  of 
chapter  IX  have  been  further  supported  by  the  two-dimensional  test  cases 
of  this  chapter.  In  respect  to  xy-geometry,  specific  conclusions  are: 

1  -  The  hybrid  Gauss-Chr istof fel  (fixed)  polar  /  equal  weight 
composite  (discrete  element)  azimuthal  angular  quadrature  is  most 
efficient  and  accurate. 

2  -  The  minimum  number  of  Gauss-Chr  istof  fel  latitudes  should  be 
two,  since  with  only  one  latitude,  the  quadrature  is  not  of  high  enough 
order  to  meet  the  diffusion  limit.  Numerical  performance  indicates  the 
importance  of  this  consideration. 

3  -  The  coupling  of  angular  and  spatial  quadrature  by  the  mechanism 
of  "steered"  element  fluxes  is  highly  successful  in  ameliorating  ray 
effects,  provided  at  least  three  (azimuthal)  elements  are  used. 

4  -  The  three-point  Gauss-Legendre  rule  was  the  most  effective 
element  quadrature  tested. 

5  -  As  in  chapter  IX,  the  expectation  that  "the  harder  the  problem, 
the  more  advantageous  the  discrete  elements  method"  is  supported  by  the 
evidence  of  these  two  tests. 

6  -  The  discrete  elements  schemes  showed  better  accuracy  and  more 
consistent  convergence  toward  the  benchmark  solution  then  the  discrete 
ordinates  schemes  tested. 

7  -  An  evaluation  of  the  full  accuracy  and  cost-effectiveness  of 
the  method  will  require  the  use  of  a  spatial  quadrature  of  higher  order, 
compatible  with  the  higher  order  of  the  discrete  elements  angular 
quadrature  (in  order  to  take  full  advantage  of  the  angular-spatial 


coupling) . 


XI.  Conclusions  and  Recommendations 


The  objective  of  this  research  has  been  to  develop  and  demonstrate 
"proof  of  concept"  of  the  discrete  elements  method  of  numerical  neutral 
particle  transport.  This  chapter  summarizes  the  conclusions  which  were 
drawn  in  the  previous  chapters  and  presents  recommendations  for  use  of 
the  discrete  elements  method  and  for  further  study. 


A.  Conclusions 

1.  Sound  Theoretical  Basis 

The  discrete  elements  method  has  a  sound  basis  in  numerical 
analysis  and  transport  theory.  It  is  not  an  ad  hoc  fixup  to  discrete 
ordinates. 

2.  Practical  for  Computer  Implementation 

The  discrete  elements  algorithm  retains  the  essential 

structure  of  the  discrete  ordinates  method.  Its  requirements  for 

computer  storage  and  run  time  are  comparable  to  S^;  in  the  absence  of 

vacuum  boundaries,  L__  uses  22%  less  storage  and  run  time  than  S,_. 
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The  method  could  be  incorporated  within  existing  production  codes. 

3.  Amelioration  of  Ray  Effects 

The  method  is  effective  in  ameliorating  ray  effects,  but  is 
more  sensitive  to  quasi-ray  effects  than  discrete  ordinates.  In  this 
regard,  a  low-order  L„  method  behaves  like  a  high  order  S„  method.  Use 
of  high-order  spatial  quadrature  should  ameliorate  the  quasi-ray  effect. 


4.  Accurac 
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The  method  is  consistently  more  accurate  than  S  with  the  same 
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quadrature  set.  For  example,  G3-L._  ,  is  more  accurate  than 
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SC 


S„„  For  the  vacuum  duct  problem,  L„  was  the  most  consistently 
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accurate  and  convergent  (as  angular  or  spatial  mesh  is  refined)  method. 
5.  Cost  Effectiveness 


The  discrete  elements  method  is  the  most  cost-effective 
alternative  for  some  problems.  This  was  demonstrated  in  slab  geometry, 
but  could  not  be  fully  demonstrated  in  xy-geometry  since  the  spatial 
quadrature  dominated  the  remaining  errors  in  the  vacuum  duct  solutions. 
The  discrete  elements  method  was  more  cost-effective  them  conventional 


S  for  the  xy-geometry  duct  problem.  Both  theory  and  testing  supported 
N 


the  conjecture  that  the  more  difficult  the  problem,  the  more 


advantageous  the  discrete  elements  method  would  be,  as  compared  to  S  . 
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6.  Element  Quadrature  Rules 

Gauss-Legendre  three-point  quadrature  was  consistently  the 
most  effective  element  quadrature  rule  tested.  It  was  more  accurate  and 
cost-effective  than  Simpson's  rule,  and  as  accurate  as  Newton-Cotes 
five-point  rule,  but  less  expensive. 

7.  Angular  Quadrature  for  XY-Geometry 

The  hybrid  Gauss-Chr istof fel  (fixed)  polar  /  equal  weight 


(steered)  azimuthal  quadrature,  L  ,  was  consistently  more  accurate 
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and  less  expensive  than  the  equal  weight  (steered)  polar  /  equal  weight 


(steered)  azimuthal  quadrature,  L 


K,L* 


8.  Spatial  Quadrature 

One  approach  considered  was  to  use  the  least  expensive  spatial 
quadrature  available,  the  step  method,  for  the  auxiliary  fluxes  (used  to 


XI-2 


O  '.'  v* 


-J 


steer  the  streaming  directions  of  the  element  "main"  flux)  with  the 
consideration  that  any  accuracy  of  steering  would  be  better  than  none. 


i.e.,  S  .  and  that  the  cheap  auxiliary  calculations  could  then  give  the 
N 

best  cost-effectiveness.  This  approach  proved  ineffective. 


The  alternative  approach  was  to  use  the  highest  order,  most 
accurate  spatial  quadrature  for  both  the  auxiliary  and  main  fluxes  so 
that  with  really  precise  steering,  the  method  would  produce  sufficient 
accuracy  to  pay  for  the  expense  of  the  high  order  auxiliary 
calculations.  This  proved  to  be  the  case.  In  one-dimensional  geometry, 
the  linear  characteristic  method  produced  the  most  accuracy  and  was  more 
cost-effective  than  Gauss-Legendre  discrete  ordinates.  In  two- 
dimensional  geometry,  the  most  effective  method  tested  was  discrete 


elements  with  step  characteristic  spatial  quadrature. 

An  analysis  of  the  quasi-ray  effect  errors  in  the  L 
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solution  to  the  vacuum  duct  problem  indicated  that  the  use  of  linear 
characteristic  spatial  quadrature  for  both  the  main  and  auxiliary  fluxes 
would  reduce  or  eliminate  those  errors. 


B.  Recommended  General  Purpose  Discrete  Elements  Schemes 


From  the  experience  gained  in  testing  these  methods,  the  following 
quadratures  are  recommended  for  general  purpose  application: 


1-D  (Slab) :  L . 


2-D  (xy)  :  L, 


No  three-dimensional  codes  were  used,  but  based  on  its  two¬ 


dimensional  performance  and  on  the  expectation  that  composite  quadrature 
with  angle-space  coupling  is  effective  for  ill-behaved  distributions, 

L  or  L  is  recommended  for  use  in  3-D  problems.  These  methods 
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should  all  use  Gauss-Legendre  three-point  element  quadrature  and  the 


highest  order  spatial  quadrature  available. 


C.  Recommended  Topics  for  Further  Research 


This  initial  study  has  developed  the  discrete  elements  method  and 
demonstrated  some  of  its  potential.  Of  necessity,  many  areas  of  interest 
were  not  pursued.  The  following  topics  are  recommended  as  being  of  most 
value  in  developing  the  discrete  elements  method  as  a  useful  tool. 


1.  Evaluation  with  Higher  Order  Spatial  Quadratures 


A  cost-effectiveness  study  of  discrete  elements  in  two¬ 


dimensional  geometry  with  linear  characteristic  or  linear  nodal  spatial 


quadrature  and  a  range  of  test  problems  could  further  explore  the  value 


of  the  method. 


2.  Extension  to  Curvilinear  Coordinates 

Characteristic  quadratures  have  not  been  applied  in 
curvilinear  coordinates.  Step  and  diamond  difference  methods  have  been, 
but  the  former  is  insufficiently  accurate  and  the  latter  is 
insufficiently  smooth  for  use  in  discrete  elements.  Extension  of  the 
method  to  curvilinear  coordinates  could  be  possible  using  the  linear 
discontinuous  spatial  quadrature. 

3.  Extension  to  Anisotropic  Scatter 

Since  the  zeroth  and  first  angular  moments  of  the  flux  are 
computed  and  used  in  the  discrete  ordinates  method,  including  linearly 
anisotropic  scatter  should  be  straightforward.  Higher  order  anisotropy 
could  be  included  by  approximating  the  higher  moments  by  numerical 
quadrature  over  each  element  as  is  done  for  the  first  moment,  or 
alternatively,  the  full  collection  of  auxiliary  fluxes  could  be  used  as 
a  single  quadrature  set  to  approximate  the  higher  angular  moments. 
Research  is  needed  to  find  the  most  accurate  and  least  expensive  scheme. 


The  discrete  elements  method  developed  here  couples  the 
angular  distribution  to  the  spatial  quadrature  only  through  the  mean 
streaming  direction.  In  a  sense,  this  is  the  first  order  member  of  a 
family  of  discrete  element  methods,  with  discrete  ordinates  as  the 
zeroth  order  member.  Higher  order  schemes  could  possibly  be  developed  by 
treating  the  flux  distribution  over  each  element  of  angle  as  a  known 
polynomial  (from  the  auxiliary  fluxes)  and  integrating  that  distribution 
across  the  space  cell  analytically,  as  an  element  characteristic  spatial 
quadrature  rather  than  the  (steered)  ordinate  characteristic  spatial 
quadratures  used  here.  Using  the  symbolic  tensor  algebra  capabilities  of 
the  MACSYMA  system,  such  methods  could  conceivably  be  extended  to 
curvilinear  coordinates.  Such  higher  order  schemes  would  abandon  some  of 
the  algorithmic  simplicity  of  Sv,  and  L  .  but  could  allow  the  use  of 
coarse  spatial  grids  and  still  provide  acceptable  accuracy.  Such  a 
hybrid  of  discrete  ordinates,  finite  elements,  and  space-angle  synthesis 
might  prove  highly  effective. 
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Abstract 


**  A  new  "discrete  elements"  (Ljj)  transport  method  is  derived  and  com¬ 
pared  to  the  discrete  ordinates  SN  method,  theoretically  and  by  numerical 
experimentation.  The  discrete  elements  method  is  more  accurate  than  dis¬ 
crete  ordinates  and  strongly  ameliorates  ray  effects  for  the  practical 
problems  studied.  The  discrete  elements  method  is  shown  to  be  more  cost 
effective,  in  terms  of  execution  time  with  comparable  storage  to  attain 
the  same  accuracy,  for  a  one-dimensional  test  case  using  linear  char¬ 
acteristic  spatial  quadrature.  In  a  two-dimensional  test  case,  a  vacuum 
duct  in  a  shield,  Ljj  is  more  consistently  convergent  toward  a  Monte  Carlo 
benchmark  solution  chan  S$j,  using  step  characteristic  spatial  quadrature. 

An  analysis  of  the  interaction  of  angular  and  spatial  quadrature  in  xy- 
geometry  indicates  the  desirability  of  using  linear  characteristic 
spatial  quadrature  with  the  Ljg  method. 

The  discrete  elements  method  is  based  on  discretizing  the  Boltzmann 
equation  over  a  set  of  elements  of  angle.  The  zeroth  and  first  angular 
moments  of  the  directional  flux,  over  each  element*,  are  estimated  by 
numerical  quadrature  and  yield  a  flux-weighted  average  streaming  direction 
for  the  element.  (Data  for  this  estimation  are  fluxes  in  fixed  directions 
calculated  as  in  S^) .  The  spatial  quadrature  then  propagates  the  element 
flux  in  this  "steered"  direction.  Since  the  quadrature  directions  are  not 
fixed,  but  are  coupled  to  the  fluxes,  the  method  strongly  ameliorates  ray 
effect.  This  is  verified  using  the  square-in-a-square  test  case  originated 
by  Lathrop.  A  variety  of  spatial,  angular,  and  element  quadrature  schemes 
are  evaluated  for  both  Ly  and  S^.  The  best  discrete  elements  method  uses 
a  hybrid  of  Gauss-Christoffel  polar  and  composite  3-point  Gauss-Legendre 
azimuthal  quadrature. 
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